Flux_parietal_Kommajosyula.cpp

/****************************************************************************
* Copyright (c) 2021, 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 <Flux_parietal_Kommajosyula.h>
#include <Hibiki_Ishii_nucleation_site_density.h>
#include <Flux_parietal_adaptatif.h>
#include <Loi_paroi_adaptative.h>
#include <Pb_Multiphase.h>
 
#include <Domaine_VF.h>
#include <TRUSTTrav.h>
#include <Milieu_composite.h>
#include <Saturation_base.h>
#include <cmath>
 
Implemente_instanciable(Flux_parietal_Kommajosyula, "Flux_parietal_Kommajosyula", Flux_parietal_base);
 
Sortie& Flux_parietal_Kommajosyula::printOn(Sortie& os) const
{
  return Flux_parietal_base::printOn(os);
}
 
Entree& Flux_parietal_Kommajosyula::readOn(Entree& is)
{
  const Pb_Multiphase& pbm = ref_cast(Pb_Multiphase, pb_.valeur());
  Correlation_base::typer_lire_correlation(correlation_monophasique_, pbm, "Flux_parietal", is);
  Cout << que_suis_je() << " : single-phase wall heat flux is " << correlation_monophasique_->que_suis_je() << finl;
 
  Param param(que_suis_je());
  param.ajouter("contact_angle_deg", &theta_, Param::REQUIRED);
  param.ajouter("molar_mass", &molar_mass_, Param::REQUIRED);
  param.lire_avec_accolades(is);
 
  if (!sub_type(Milieu_composite, pb_->milieu()))
    Process::exit("Flux_parietal_Kommajosyula::readOn: the medium must be composite!");
 
  if (!pbm.nom_phase(0).debute_par("liquide"))
    Process::exit("Flux_parietal_Kommajosyula::readOn: the first phase must be liquid!");
 
  const Milieu_composite& milc = ref_cast(Milieu_composite, pb_->milieu());
  int n_g = -1;
  for (int k = 1; k < pbm.nb_phases() ; k++)
    if (milc.has_saturation(0, k))
      n_g += 1;
 
  if (n_g > 0)
    Process::exit("Flux_parietal_Kommajosyula::readOn: there can only be one evaporating phase for the carrying liquid for now! Please feel free to update the code if you need.");
 
  return is;
}
 
void Flux_parietal_Kommajosyula::completer()
{
  correlation_monophasique_->completer();
}
 
void Flux_parietal_Kommajosyula::qp(const input_t& in, output_t& out) const
{
  // Set everything to zero (or one)
  if (out.qpk)     (*out.qpk)     = 0.;
  if (out.da_qpk)  (*out.da_qpk)  = 0.;
  if (out.dp_qpk)  (*out.dp_qpk)  = 0.;
  if (out.dv_qpk)  (*out.dv_qpk)  = 0.;
  if (out.dTf_qpk) (*out.dTf_qpk) = 0.;
  if (out.dTp_qpk) (*out.dTp_qpk) = 0.;
  if (out.qpi)     (*out.qpi)     = 0.;
  if (out.da_qpi)  (*out.da_qpi)  = 0.;
  if (out.dp_qpi)  (*out.dp_qpi)  = 0.;
  if (out.dv_qpi)  (*out.dv_qpi)  = 0.;
  if (out.dTf_qpi) (*out.dTf_qpi) = 0.;
  if (out.dTp_qpi) (*out.dTp_qpi) = 0.;
  if (out.nonlinear) (*out.nonlinear) = 1;
 
  // Single phase, no interfacial flux
  ref_cast(Flux_parietal_base, correlation_monophasique_.valeur()).qp(in, out);
 
  // The liquid phase is necessarily phase 0 because the single-phase correlation only fills phase 0
  int n_l = 0;
  const Milieu_composite& milc = ref_cast(Milieu_composite, pb_->milieu());
 
  for (int k = 0; k < in.N; k++)
    if (n_l != k)
      if (milc.has_saturation(n_l, k))
        {
          const int ind_sat = k < n_l
                              ? (k * (in.N - 1) - (k - 1) * k / 2) + (n_l - k - 1)
                              : (n_l * (in.N - 1) - (n_l - 1) * n_l / 2) + (k - n_l - 1);
 
          const double Delta_T_sup = in.Tp - in.Tsat[ind_sat]; // Wall superheat
 
          if (Delta_T_sup > 0) // No wall superheat => no nucleation => single phase heat transfer only
            {
              const double Delta_T_sub = std::max(in.Tsat[ind_sat] - in.T[n_l], 1.e-8); // Subcooling (non-negative for numerical reasons)
 
              double u_bulk = 0;
              if (sub_type(Flux_parietal_adaptatif, correlation_monophasique_.valeur()))
                {
                  const Loi_paroi_adaptative& corr_loi_paroi = ref_cast(Loi_paroi_adaptative, ref_cast(Pb_Multiphase, pb_.valeur()).get_correlation("Loi_paroi"));
                  const double u_tau = corr_loi_paroi.get_utau(in.f);
                  u_bulk = 20. * u_tau; // WARNING: Big approximation...
                }
              else
                Cerr << "Flux_parietal_Kommajosyula::qp isn't adapted to a single-phase " << correlation_monophasique_->que_suis_je() << " heat flux for now !\n", Process::exit();
 
              // Single phase heat transfer coefficient
              const double h_fc = (*out.dTp_qpk)(n_l);
 
              // Compute nucleation parameters (diameters, frequencies, site density)
              const NucleationParams nuc = compute_nucleation_params(in, n_l, k, ind_sat, Delta_T_sup, Delta_T_sub, u_bulk);
 
              // Thermal boundary layer reformation time (page 61 Kommajosyula PhD; error in equation 4.5 of the manuscript)
              const double t_star = in.lambda[n_l] * in.rho[n_l] * in.Cp[n_l] / (h_fc * h_fc * M_PI);
 
              // Compute sliding area and surface fraction
              const SlidingAreaParams sa = compute_sliding_area(nuc, h_fc, t_star);
 
              DoubleTrav qpk_loc, da_qpk_loc, dp_qpk_loc, dv_qpk_loc, dTf_qpk_loc, dTp_qpk_loc;
 
              if (out.qpk)     qpk_loc     = *out.qpk;
              if (out.da_qpk)  da_qpk_loc  = *out.da_qpk;
              if (out.dp_qpk)  dp_qpk_loc  = *out.dp_qpk;
              if (out.dv_qpk)  dv_qpk_loc  = *out.dv_qpk;
              if (out.dTf_qpk) dTf_qpk_loc = *out.dTf_qpk;
              if (out.dTp_qpk) dTp_qpk_loc = *out.dTp_qpk;
 
              // Correct the single phase heat flux
              if (out.qpk)     (*out.qpk)    *= (1 - sa.S_sl);
              if (out.da_qpk)  (*out.da_qpk) *= (1 - sa.S_sl);
              if (out.dp_qpk)  (*out.dp_qpk) *= (1 - sa.S_sl);
              if (out.dv_qpk)  (*out.dv_qpk) *= (1 - sa.S_sl);
              if (out.dTf_qpk)
                {
                  (*out.dTf_qpk)          *= (1 - sa.S_sl);
                  (*out.dTf_qpk)(n_l,n_l) += -sa.dTl_S_sl * qpk_loc(n_l);
                  (*out.dTf_qpk)(n_l, k ) += -sa.dTk_S_sl * qpk_loc( k );
                }
              if (out.dTp_qpk)
                {
                  (*out.dTp_qpk)      *= (1 - sa.S_sl);
                  (*out.dTp_qpk)(n_l) += -sa.dTp_S_sl * qpk_loc(n_l);
                }
 
              // Bubble sliding
              if (out.qpk)     (*out.qpk)(n_l)           += sa.S_sl * 2. * h_fc * (in.Tp - in.T[n_l]);
              if (out.dTf_qpk) (*out.dTf_qpk)(n_l, n_l)  += -sa.S_sl * 2. * h_fc + sa.dTl_S_sl * 2. * h_fc * (in.Tp - in.T[n_l]);
              if (out.dTf_qpk) (*out.dTf_qpk)(n_l, k)    += sa.dTk_S_sl * 2. * h_fc * (in.Tp - in.T[n_l]);
              if (out.dTp_qpk) (*out.dTp_qpk)( k )       += sa.S_sl * 2. * h_fc + sa.dTp_S_sl * 2. * h_fc * (in.Tp - in.T[n_l]);
 
              // Evaporation
              if (out.qpi)         (*out.qpi)(n_l, k) = 1./6. * M_PI * in.rho[k] * in.Lvap[ind_sat] *                    std::pow(nuc.D_lo, 3.) *         nuc.f_dep *     sa.N_active;
              if (out.dTp_qpi) (*out.dTp_qpi)(n_l, k) = 1./6. * M_PI * in.rho[k] * in.Lvap[ind_sat] * (3. * nuc.dTp_D_lo * std::pow(nuc.D_lo, 2.) *         nuc.f_dep *     sa.N_active
                                                                                                         +                std::pow(nuc.D_lo, 3.) * nuc.dTp_f_dep *     sa.N_active
                                                                                                         +                std::pow(nuc.D_lo, 3.) *     nuc.f_dep * sa.dTp_N_active);
              if (out.dTf_qpi) (*out.dTf_qpi)(n_l, k, n_l) = 1./6. * M_PI * in.rho[k] * in.Lvap[ind_sat] * (3. * nuc.dTl_D_lo * std::pow(nuc.D_lo, 2.) *         nuc.f_dep *     sa.N_active
                                                                                                              +                std::pow(nuc.D_lo, 3.) * nuc.dTl_f_dep *     sa.N_active
                                                                                                              +                std::pow(nuc.D_lo, 3.) *     nuc.f_dep * sa.dTl_N_active);
              if (out.dTf_qpi) (*out.dTf_qpi)(n_l, k, k) = 1./6. * M_PI * in.rho[k] * in.Lvap[ind_sat] * (std::pow(nuc.D_lo, 3.) * nuc.f_dep * sa.dTk_N_active);
 
              if (out.d_nuc) (*out.d_nuc)(k) = nuc.D_lo;
 
            }
        }
}
 
Flux_parietal_Kommajosyula::NucleationParams
Flux_parietal_Kommajosyula::compute_nucleation_params(const input_t& in, int n_l, int k, int ind_sat,
                                                      double Delta_T_sup, double Delta_T_sub, double u_bulk) const
{
  NucleationParams nuc {};
 
  const double dTp_Delta_T_sup = 1.;
  const double dTl_Delta_T_sub = -1.;
  const double Ja_sup = in.rho[n_l] * in.Cp[n_l] * Delta_T_sup / (in.rho[k] * in.Lvap[ind_sat]);
  const double Ja_sub = in.rho[n_l] * in.Cp[n_l] * Delta_T_sub / (in.rho[k] * in.Lvap[ind_sat]);
  const double dTp_Ja_sup = in.rho[n_l] * in.Cp[n_l] / (in.rho[k] * in.Lvap[ind_sat]);
  const double dTl_Ja_sub = -in.rho[n_l] * in.Cp[n_l] / (in.rho[k] * in.Lvap[ind_sat]);
 
  // Nucleation site density (Hibiki Ishii 2003) — T_ref = Tp (wall temperature)
  nuc.N_sites     =     Hibiki_Ishii_site_density(in.rho[k], in.rho[n_l], in.Tp, in.p, in.Lvap[ind_sat], in.Tsat[ind_sat], in.Sigma[ind_sat], theta_, molar_mass_);
  nuc.dTp_N_sites = dT_ref_Hibiki_Ishii_site_density(in.rho[k], in.rho[n_l], in.Tp, in.p, in.Lvap[ind_sat], in.Tsat[ind_sat], in.Sigma[ind_sat], theta_, molar_mass_);
  nuc.dTl_N_sites = 0.;
  nuc.dTk_N_sites = 0.;
 
  // Departure diameter (page 47 Kommajosyula PhD)
  const double rho_ratio_pow = std::pow((in.rho[n_l] - in.rho[k]) / in.rho[k], 0.27);
  const double u_bulk_pow = std::pow(u_bulk, -0.26);
  nuc.D_d =
    18.9e-6 * rho_ratio_pow * std::pow(Ja_sup, 0.75) * std::pow(1 + Ja_sub, -0.3) * u_bulk_pow;
  nuc.dTp_D_d = 18.9e-6 * rho_ratio_pow * 0.75 * dTp_Ja_sup * std::pow(Ja_sup, -0.25) *
                std::pow(1 + Ja_sub, -0.3) * u_bulk_pow;
  nuc.dTl_D_d = 18.9e-6 * rho_ratio_pow * std::pow(Ja_sup, 0.75) * -0.3 * dTl_Ja_sub *
                std::pow(1 + Ja_sub, -1.3) * u_bulk_pow;
 
  // Liftoff diameter (page 47 Kommajosyula PhD)
  nuc.D_lo     = 1.2 * nuc.D_d;
  nuc.dTp_D_lo = 1.2 * nuc.dTp_D_d;
  nuc.dTl_D_lo = 1.2 * nuc.dTl_D_d;
 
  // Bubble growth constant (page 42 Kommajosyula PhD) excluding microlayer contribution
  const double chi     = 1.55 - 0.05 * Delta_T_sub / Delta_T_sup;
  const double dTp_chi = 0.05 * Delta_T_sub * dTp_Delta_T_sup / std::pow(Delta_T_sup, 2);
  const double dTl_chi = -0.05 * dTl_Delta_T_sub / Delta_T_sup;
  const double thermal_diff = in.lambda[n_l] / (in.rho[n_l] * in.Cp[n_l]);
  const double growth_factor = 2. * std::sqrt(3. / M_PI) * thermal_diff;
 
  const double K     = chi     * growth_factor * Ja_sup;
  const double dTp_K = dTp_chi * growth_factor * Ja_sup + chi * growth_factor * dTp_Ja_sup;
  const double dTl_K = dTl_chi * growth_factor * Ja_sup;
 
  // Bubble growth time (page 44 Kommajosyula PhD)
  nuc.t_g     = nuc.D_d * nuc.D_d / (16. * K * K);
  nuc.dTp_t_g = 2 * nuc.dTp_D_d * nuc.D_d / (16. * K * K)
                + nuc.D_d * nuc.D_d * -2. * dTp_K / (16. * K * K * K);
  nuc.dTl_t_g = 2 * nuc.dTl_D_d * nuc.D_d / (16. * K * K)
                + nuc.D_d * nuc.D_d * -2. * dTl_K / (16. * K * K * K);
 
  // Bubble wait time (page 56 Kommajosyula PhD)
  const double t_w     = 0.0061 * std::pow(Ja_sub, 0.63) / Delta_T_sup;
  const double dTp_t_w = 0.0061 * std::pow(Ja_sub, 0.63) * -dTp_Delta_T_sup / std::pow(Delta_T_sup, 2);
  const double dTl_t_w = 0.0061 * .63 * dTl_Ja_sub * std::pow(Ja_sub, -.37) / Delta_T_sup;
 
  // Bubble departure frequency (page 51 Kommajosyula PhD)
  nuc.f_dep     =  1 / (nuc.t_g + t_w);
  nuc.dTp_f_dep = -(nuc.dTp_t_g + dTp_t_w) / std::pow(nuc.t_g + t_w, 2);
  nuc.dTl_f_dep = -(nuc.dTl_t_g + dTl_t_w) / std::pow(nuc.t_g + t_w, 2);
 
  return nuc;
}
 
Flux_parietal_Kommajosyula::SlidingAreaParams
Flux_parietal_Kommajosyula::compute_sliding_area(const NucleationParams& nuc, double h_fc, double t_star) const
{
  SlidingAreaParams sa {};
 
  // Active nucleation site density (page 66 Kommajosyula PhD)
  const double N_0     = nuc.f_dep     * nuc.t_g     * M_PI * nuc.D_d * nuc.D_d / 4.;
  const double dTp_N_0 = nuc.dTp_f_dep * nuc.t_g     * M_PI * nuc.D_d * nuc.D_d / 4.
                         + nuc.f_dep   * nuc.dTp_t_g * M_PI * nuc.D_d * nuc.D_d / 4.
                         + nuc.f_dep   * nuc.t_g     * M_PI * 2. * nuc.dTp_D_d * nuc.D_d / 4.;
  const double dTl_N_0 = nuc.dTl_f_dep * nuc.t_g     * M_PI * nuc.D_d * nuc.D_d / 4.
                         + nuc.f_dep   * nuc.dTl_t_g * M_PI * nuc.D_d * nuc.D_d / 4.
                         + nuc.f_dep   * nuc.t_g     * M_PI * 2. * nuc.dTl_D_d * nuc.D_d / 4.;
 
  const double W_val     =    Lambert_W_function(N_0 * nuc.N_sites);
  const double dLoc_Lambert = dx_Lambert_W_function(N_0 * nuc.N_sites);
 
  sa.N_active     = W_val / N_0;
  sa.dTp_N_active = (dTp_N_0 * nuc.N_sites + N_0 * nuc.dTp_N_sites) * dLoc_Lambert / N_0
                    + W_val * -dTp_N_0 / std::pow(N_0, 2);
  sa.dTl_N_active = (dTl_N_0 * nuc.N_sites + N_0 * nuc.dTl_N_sites) * dLoc_Lambert / N_0
                    + W_val * -dTl_N_0 / std::pow(N_0, 2);
  sa.dTk_N_active = N_0 * nuc.dTk_N_sites * dLoc_Lambert / N_0;
 
  if (sa.N_active < 0)
    Cerr << "nsite " << nuc.N_sites << "N0 " << N_0 << "\n";
 
  // Single bubble sliding area
  const double A_sl = std::pow(sa.N_active, -.5) * (nuc.D_d + nuc.D_lo) / 2.;
  const double dTp_A_sl =
    -.5 * sa.dTp_N_active * std::pow(sa.N_active, -1.5) * (nuc.D_d + nuc.D_lo) / 2. +
    std::pow(sa.N_active, -.5) * (nuc.dTp_D_d + nuc.dTp_D_lo) / 2.;
  const double dTl_A_sl =
    -.5 * sa.dTl_N_active * std::pow(sa.N_active, -1.5) * (nuc.D_d + nuc.D_lo) / 2. +
    std::pow(sa.N_active, -.5) * (nuc.dTl_D_d + nuc.dTl_D_lo) / 2.;
  const double dTk_A_sl =
    -.5 * sa.dTk_N_active * std::pow(sa.N_active, -1.5) * (nuc.D_d + nuc.D_lo) / 2.;
 
  // Percentage of the surface occupied by sliding bubbles
  const double product = A_sl * sa.N_active * nuc.f_dep * t_star;
  const bool saturated = (1. < product);
  sa.S_sl     = std::min(1., product);
  sa.dTp_S_sl = saturated ? 0. :
                dTp_A_sl    * sa.N_active     * nuc.f_dep     * t_star
                + A_sl      * sa.dTp_N_active * nuc.f_dep     * t_star
                + A_sl      * sa.N_active     * nuc.dTp_f_dep * t_star;
  sa.dTl_S_sl = saturated ? 0. :
                dTl_A_sl    * sa.N_active     * nuc.f_dep     * t_star
                + A_sl      * sa.dTl_N_active * nuc.f_dep     * t_star
                + A_sl      * sa.N_active     * nuc.dTl_f_dep * t_star;
  sa.dTk_S_sl = saturated ? 0. :
                dTk_A_sl    * sa.N_active     * nuc.f_dep     * t_star
                + A_sl      * sa.dTk_N_active * nuc.f_dep     * t_star;
 
  return sa;
}
 
double Flux_parietal_Kommajosyula::Lambert_W_function(double x) const
{
  if (x < 1/M_E)
    return x;
  else if (x > M_E)
    return std::log(x) - std::log(std::log(x));
  else
    return 0.2689*x + 0.2690;
}
 
double Flux_parietal_Kommajosyula::dx_Lambert_W_function(double x) const
{
  if (x < 1/M_E)
    return 1.;
  else if (x > M_E)
    return 1./x - 1./(x*std::log(x));
  else
    return 0.2689;
}