Loi_Etat_Multi_GP_WC.cpp

/****************************************************************************
* Copyright (c) 2026, CEA
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
* 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
* 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
* OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
*****************************************************************************/
 
#include <Fluide_Weakly_Compressible.h>
#include <Loi_Etat_Multi_GP_WC.h>
#include <Champ_Uniforme.h>
 
#include <Probleme_base.h>
#include <Domaine_VF.h>
#include <Param.h>
#include <Debog.h>
 
Implemente_instanciable_sans_constructeur(Loi_Etat_Multi_GP_WC,"Loi_Etat_Multi_Gaz_Parfait_WC",Loi_Etat_Multi_GP_base);
// XD multi_gaz_parfait_WC loi_etat_gaz_parfait_base multi_gaz_parfait_WC INHERITS_BRACE Class for perfect gas
// XD_CONT multi-species mixtures state law used with a weakly-compressible fluid.
 
Loi_Etat_Multi_GP_WC::Loi_Etat_Multi_GP_WC() : num_espece_(-1) { }
 
Sortie& Loi_Etat_Multi_GP_WC::printOn(Sortie& os) const
{
  os <<que_suis_je()<< finl;
  return os;
}
 
Entree& Loi_Etat_Multi_GP_WC::readOn(Entree& is)
{
  Param param(que_suis_je());
  param.ajouter("Species_number",&num_espece_,Param::REQUIRED); // XD_ADD_P entier
  // XD_CONT Number of species you are considering in your problem.
  param.ajouter("Diffusion_coeff",&diffusion_coeff_,Param::REQUIRED); // XD_ADD_P field_base
  // XD_CONT Diffusion coefficient of each species, defined with a Champ_uniforme of dimension equals to the
  // XD_CONT species_number.
  param.ajouter("Molar_mass",&molar_mass_,Param::REQUIRED); // XD_ADD_P field_base
  // XD_CONT Molar mass of each species, defined with a Champ_uniforme of dimension equals to the species_number.
  param.ajouter("Mu",&mu_,Param::REQUIRED); // XD_ADD_P field_base
  // XD_CONT Dynamic viscosity of each species, defined with a Champ_uniforme of dimension equals to the species_number.
  param.ajouter("Cp",&cp_,Param::REQUIRED);// XD_ADD_P field_base
  // XD_CONT Specific heat at constant pressure of the gas Cp, defined with a Champ_uniforme of dimension equals to the
  // XD_CONT species_number..
  param.ajouter("Prandtl",&Pr_,Param::REQUIRED); // XD_ADD_P double
  // XD_CONT Prandtl number of the gas Pr=mu*Cp/lambda.
  param.lire_avec_accolades_depuis(is);
 
  // XXX
  if (diffusion_coeff_->valeurs().line_size() != num_espece_ ||
      molar_mass_->valeurs().line_size() != num_espece_ ||
      mu_->valeurs().line_size() != num_espece_  ||
      cp_->valeurs().line_size() != num_espece_ )
    {
      Cerr << "Error while reading the EOS Multi_Gaz_Parfait_WC : This is not allowed !" << finl;
      Cerr << "Verify that the dimension of the fields Diffusion_coeff, Molar_mass, Mu and Cp is = " << num_espece_ << finl;
      Process::exit();
    }
 
  return is;
}
 
void Loi_Etat_Multi_GP_WC::initialiser_inco_ch()
{
  /*
   * XXX : In this problem we solve n-1 equations for n species
   * The mass fraction of species n is [ 1 - (y1+y2+ ... +yn-1) ]
   */
 
  if (liste_Y.size() + 1  != num_espece_ )
    {
      Cerr << "Error in your data file : This is not allowed !" << finl;
      Cerr << "You should define " << num_espece_-1 << " Species equations and not " << liste_Y.size() << finl;
      Process::exit();
    }
 
  Fluide_Weakly_Compressible& FWC = ref_cast(Fluide_Weakly_Compressible,le_fluide.valeur());
  Champ_Don_base& Yn = FWC.fraction_massique_nonresolue();
  Yn.nommer("fraction_massique_nonresolue");
  double t = le_fluide->masse_volumique().temps(); // not 0 because resuming is possible
  update_Yn_values(Yn,t);
  Loi_Etat_Multi_GP_base::initialiser_inco_ch();
}
 
void Loi_Etat_Multi_GP_WC::update_Yn_values(Champ_Don_base& Yn, double temps)
{
  DoubleTab& tab_Yn = Yn.valeurs();
  Yn.mettre_a_jour(temps);
  tab_Yn = 1.0;// re-initialize as 1
  const int size = liste_Y(0)->valeurs().size();
  assert(tab_Yn.size() == size);
  assert (liste_Y.size() + 1  == num_espece_);
 
  for (int i=0; i<liste_Y.size(); i++)
    for (int n=0; n<size; n++) tab_Yn(n,0) -= liste_Y(i)->valeurs()(n,0);
 
  double min_Yn = local_min_vect(tab_Yn);
  if (min_Yn < 0)
    {
      Cerr << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << finl;
      Cerr << "BIG PROBLEM : check your calculation because we are detecting negative mass fraction !" << finl;
      Cerr << "The results are saved to help you !" << finl;
      Cerr << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << finl;
      Fluide_Weakly_Compressible& FWC = ref_cast(Fluide_Weakly_Compressible,le_fluide.valeur());
      Probleme_base& pb = FWC.get_problem();
      pb.postraiter(1);
      pb.sauver();
      Process::exit();
    }
}
 
/*! @brief Computes the molar mass of the mixture (M). M depends on the molar mass of each species and the mixture composition (Yi).
 *
 */
void Loi_Etat_Multi_GP_WC::calculer_masse_molaire(DoubleTab& tab_masse_mol_mel) const
{
  assert (sub_type(Champ_Uniforme,masse_molaire_especes()));
 
  const int size = tab_masse_mol_mel.size();
  assert (liste_Y(0)->valeurs().size()==size);
  ArrOfDouble inv_M(size);
 
  const Fluide_Weakly_Compressible& FWC = ref_cast(Fluide_Weakly_Compressible,le_fluide.valeur());
  const DoubleTab& Yn = FWC.fraction_massique_nonresolue().valeurs();
 
  for (int i=0; i < num_espece_; i++)
    {
      const DoubleTab& Y_i = (i == num_espece_ -1) ? Yn : liste_Y(i)->valeurs();
      const double M_i = masse_molaire_especes().valeurs()(0,i);
      for (int elem=0; elem<size; elem++) inv_M[elem] += Y_i(elem,0)/M_i;
    }
 
  for (int elem=0; elem<size; elem++) tab_masse_mol_mel(elem,0) = 1.0 / inv_M[elem];
}
 
/*! @brief Computes the Cp of the mixture. Cp depends on the Cp of each species and the mixture composition (Yi).
 *
 */
void Loi_Etat_Multi_GP_WC::calculer_tab_Cp(DoubleTab& tab_Cp) const
{
  // FIXME : Currently Cp is assumed constant for each species
  assert (sub_type(Champ_Uniforme,cp_especes()));
 
  tab_Cp = 0;
  const int size =liste_Y(0)->valeurs().size();
 
  const Fluide_Weakly_Compressible& FWC = ref_cast(Fluide_Weakly_Compressible,le_fluide.valeur());
  const DoubleTab& Yn = FWC.fraction_massique_nonresolue().valeurs();
 
  for (int i=0; i<num_espece_; i++)
    {
      // TODO : FIXME : check whether Yn.valeurs() is up to date
      const DoubleTab& Y_i = (i == num_espece_ -1) ? Yn : liste_Y(i)->valeurs();
      const double cp_i = cp_especes().valeurs()(0,i);
      for (int elem=0; elem<size; elem++) tab_Cp(elem,0) += Y_i(elem,0)*cp_i;
    }
}
 
/*! @brief Recomputes the density.
 *
 */
void Loi_Etat_Multi_GP_WC::calculer_masse_volumique()
{
  const DoubleTab& tab_ICh = le_fluide->inco_chaleur().valeurs();
  DoubleTab& tab_rho = le_fluide->masse_volumique().valeurs();
 
  Fluide_Weakly_Compressible& FWC = ref_cast(Fluide_Weakly_Compressible,le_fluide.valeur());
  Champ_Don_base& Yn = FWC.fraction_massique_nonresolue();
  double temps = le_fluide->masse_volumique().temps();
  update_Yn_values(Yn,temps); // XXX : check whether the call is in the right place ...
  Loi_Etat_Multi_GP_base::calculer_masse_molaire();
 
  assert (tab_rho.line_size() == 1);
  const int n = tab_rho.size();
  const DoubleTab& pres = FWC.pression_th_tab();
 
  for (int som=0 ; som<n ; som++)
    {
      double r = R_GAS / masse_mol_mel(som);
      tab_rho_np1(som) = calculer_masse_volumique(pres(som,0),tab_ICh(som,0),r);
      tab_rho(som,0) = 0.5 * ( tab_rho_n(som) + tab_rho_np1(som) );
    }
 
  tab_rho.echange_espace_virtuel();
  tab_rho_np1.echange_espace_virtuel();
  le_fluide->calculer_rho_face(tab_rho_np1);
  Debog::verifier("Loi_Etat_Multi_GP_WC::calculer_masse_volumique, tab_rho_np1",tab_rho_np1);
  Debog::verifier("Loi_Etat_Multi_GP_WC::calculer_masse_volumique, tab_rho",tab_rho);
}
 
double Loi_Etat_Multi_GP_WC::calculer_masse_volumique(double P, double T, double r) const
{
  return Loi_Etat_Multi_GP_base::calculer_masse_volumique(P,T,r);
}
 
double Loi_Etat_Multi_GP_WC::calculer_masse_volumique(double P, double T) const
{
  return Loi_Etat_Multi_GP_base::calculer_masse_volumique(P,T);
}
 
/*! @brief Computes the reference dynamic viscosity (depends on mass fractions Yi).
 *
 * With Wilke formulation: https://aip.scitation.org/doi/pdf/10.1063/1.1747673
 * See also for mass fractions : https://doi.org/10.1016/j.ijheatmasstransfer.2020.120470
 *
 *
 * Mu =
 *
 *       N
 *     ______
 *     \     `
 *     \          Mu_i*Y_i
 *      \     ----------------
 *       \      N
 *        \    __
 *        /    \ `
 *       /      )   Phi_ij*Y_j
 *      /      /_,
 *     /      j = 1
 *     /_____,
 *     i = 1
 *
 *
 * where Phi_ij =
 *
 *                                   2
 *         /     0.25     ______    \
 *         |/M_j\        / Mu_i     |
 *     M_i*||---|    *  /  ----  + 1|
 *         \\M_i/     \/   Mu_j     /
 *     -------------------------------
 *                   ___________
 *                  / 8*M_i
 *           M_j*  /  ----- + 8
 *               \/    M_j
 *
 *
 */
void Loi_Etat_Multi_GP_WC::calculer_mu_wilke()
{
  assert (liste_Y.size() + 1 == num_espece_);
  const int size = liste_Y(0)->valeurs().size();
  DoubleTrav phi(size), mu(size);
  mu = 0.;
  phi = 0.;
 
  Fluide_Weakly_Compressible& FWC = ref_cast(Fluide_Weakly_Compressible,le_fluide.valeur());
  DoubleTab& Yn = FWC.fraction_massique_nonresolue().valeurs();
 
  for (int i=0; i < num_espece_ ; i++)
    {
      phi = 0.;
      const double M_i = masse_molaire_especes().valeurs()(0,i);
      const double mu_i = visc_dynamique_especes().valeurs()(0,i);
 
      for (int j=0; j < num_espece_ ; j++)
        if (j != i) // otherwise phi_ii = 1
          {
            const double M_j = masse_molaire_especes().valeurs()(0,j);
            const double mu_j = visc_dynamique_especes().valeurs()(0,j);
 
            double a = 1. + sqrt( mu_i / mu_j ) * pow( M_j / M_i , 0.25);
            double b = sqrt( 8. * ( 1. + ( M_i / M_j )));
            double phi_ij = ( M_i / M_j ) * a * a / b;
            // TODO : FIXME : check whether Yn.valeurs() is up to date
            const DoubleTab& y_j = (j == num_espece_ -1) ? Yn : liste_Y(j)->valeurs();
            // node is elem (VDF) or face (VEF)
            for (int node=0; node<y_j.size(); node++) phi(node) += y_j(node,0) * phi_ij;
          }
      // We add the mass fraction when i = j
      const DoubleTab& y_i = (i == num_espece_ -1) ? Yn : liste_Y(i)->valeurs();
      for (int node=0; node<y_i.size(); node++) mu(node) += mu_i * y_i(node,0) / ( y_i(node,0) + phi(node) );
    }
 
  calculer_tab_mu(mu, size);
}
 
/*! @brief Computes the dynamic viscosity divided by the Schmidt number = rho * D.
 *
 */
void Loi_Etat_Multi_GP_WC::calculer_mu_sur_Sc()
{
  /*
   * ====================================================================
   * Schmidt number : Sc = nu / D = mu / ( rho * D )
   * => rho * D = mu/Sc
   * species equation : div ( rho * D * grad Y1) = div ( mu/Sc * grad Y1)
   * ====================================================================
   */
 
  Champ_Don_base& mu_sur_Sc = le_fluide->mu_sur_Schmidt();
  DoubleTab& tab_mu_sur_Sc = mu_sur_Sc.valeurs();
  const Champ_base& rho = le_fluide->masse_volumique();
  const DoubleTab& tab_rho = rho.valeurs();
  const int n = tab_mu_sur_Sc.size();
 
  // TODO : FIXME : tab_mu_sur_Sc.line_size() = 1 :( :/ :(
  // BUG : line_size should equal num_espece_ since D can vary between species
 
  if (!sub_type(Champ_Uniforme,mu_sur_Sc))
    {
      if (sub_type(Champ_Uniforme,rho))
        {
          Cerr << "We should not have a density field of type Champ_Uniforme !" << finl;
          Process::exit();
        }
      else
        for (int i=0 ; i<n ; i++)
          {
            // TODO : FIXME : D from species 1 is used here ...
            tab_mu_sur_Sc(i,0) = tab_rho(i,0) * coeff_diffusion_especes().valeurs()(0,0);
          }
    }
  else
    {
      Cerr << "We should not have a mu_sur_Sc of type Champ_Uniforme !" << finl;
      Process::exit();
    }
 
  mu_sur_Sc.changer_temps(rho.temps());
  tab_mu_sur_Sc.echange_espace_virtuel();
}
 
/*! @brief Computes the kinematic viscosity divided by the Schmidt number = D.
 *
 */
void Loi_Etat_Multi_GP_WC::calculer_nu_sur_Sc()
{
  /*
   * ====================================================================
   * Schmidt number : Sc = nu / D
   * => D = nu/Sc
   * ====================================================================
   */
 
  Champ_Don_base& nu_sur_Sc = le_fluide->nu_sur_Schmidt();
  DoubleTab& tab_nu_sur_Sc = nu_sur_Sc.valeurs();
  const int n = tab_nu_sur_Sc.size();
 
  // TODO : FIXME : tab_nu_sur_Sc.line_size() = 1 :( :/ :(
  // BUG : line_size should equal num_espece_ since D can vary between species
 
  // TODO : FIXME : D from species 1 is used here ...
  for (int i=0 ; i<n ; i++) tab_nu_sur_Sc(i,0) = coeff_diffusion_especes().valeurs()(0,0);
 
  double temps_champ = le_fluide->masse_volumique().temps();
  nu_sur_Sc.changer_temps(temps_champ);
  tab_nu_sur_Sc.echange_espace_virtuel();
}