Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "WCNSFV2PInterfaceAreaSourceSink.h"
      11             : #include "NS.h"
      12             : #include "NonlinearSystemBase.h"
      13             : #include "NavierStokesMethods.h"
      14             : 
      15             : registerMooseObject("NavierStokesApp", WCNSFV2PInterfaceAreaSourceSink);
      16             : 
      17             : InputParameters
      18         123 : WCNSFV2PInterfaceAreaSourceSink::validParams()
      19             : {
      20         123 :   InputParameters params = FVElementalKernel::validParams();
      21         123 :   params.addClassDescription(
      22             :       "Source and sink of interfacial area for two-phase flow mixture model.");
      23         246 :   params.addRequiredParam<MooseFunctorName>("u", "The velocity in the x direction.");
      24         246 :   params.addParam<MooseFunctorName>("v", "The velocity in the y direction.");
      25         246 :   params.addParam<MooseFunctorName>("w", "The velocity in the z direction.");
      26             : 
      27         246 :   params.addParam<MooseFunctorName>("L", 1.0, "The characteristic dissipation length.");
      28         123 :   params.addRequiredParam<MooseFunctorName>(NS::density, "Continuous phase density.");
      29         246 :   params.addRequiredParam<MooseFunctorName>(NS::density + std::string("_d"),
      30             :                                             "Dispersed phase density.");
      31         123 :   params.addRequiredParam<MooseFunctorName>(NS::pressure, "Continuous phase density.");
      32         246 :   params.addParam<MooseFunctorName>(
      33         246 :       "k_c", 0.0, "Mass exchange coefficients from continous to dispersed phases.");
      34         246 :   params.addParam<MooseFunctorName>("fd", 0.0, "Fraction dispersed phase.");
      35         246 :   params.addParam<Real>("fd_max", 1.0, "Maximum dispersed phase fraction.");
      36             : 
      37         246 :   params.addParam<MooseFunctorName>("sigma", 1.0, "Surface tension between phases.");
      38         246 :   params.addParam<MooseFunctorName>("particle_diameter", 1.0, "Maximum particle diameter.");
      39             : 
      40         246 :   params.addParam<Real>("cutoff_fraction",
      41         246 :                         0.1,
      42             :                         "Void fraction at which the interface area density mass transfer model is "
      43             :                         "activated. Below this fraction, spherical bubbles are assumed.");
      44             : 
      45         123 :   params.set<unsigned short>("ghost_layers") = 2;
      46             : 
      47         123 :   return params;
      48           0 : }
      49             : 
      50          66 : WCNSFV2PInterfaceAreaSourceSink::WCNSFV2PInterfaceAreaSourceSink(const InputParameters & params)
      51             :   : FVElementalKernel(params),
      52          66 :     _dim(_subproblem.mesh().dimension()),
      53         132 :     _u_var(getFunctor<ADReal>("u")),
      54         198 :     _v_var(params.isParamValid("v") ? &(getFunctor<ADReal>("v")) : nullptr),
      55          66 :     _w_var(params.isParamValid("w") ? &(getFunctor<ADReal>("w")) : nullptr),
      56         132 :     _characheristic_length(getFunctor<ADReal>("L")),
      57          66 :     _rho_mixture(getFunctor<ADReal>(NS::density)),
      58         198 :     _rho_d(getFunctor<ADReal>(NS::density + std::string("_d"))),
      59          66 :     _pressure(getFunctor<ADReal>(NS::pressure)),
      60         132 :     _mass_exchange_coefficient(getFunctor<ADReal>("k_c")),
      61         132 :     _f_d(getFunctor<ADReal>("fd")),
      62         132 :     _f_d_max(getParam<Real>("fd_max")),
      63         132 :     _sigma(getFunctor<ADReal>("sigma")),
      64         132 :     _particle_diameter(getFunctor<ADReal>("particle_diameter")),
      65         198 :     _cutoff_fraction(getParam<Real>("cutoff_fraction"))
      66             : {
      67          66 :   if (_dim >= 2 && !_v_var)
      68           0 :     paramError("v", "In two or more dimensions, the v velocity must be supplied!");
      69             : 
      70          66 :   if (_dim >= 3 && !_w_var)
      71           0 :     paramError("w", "In three or more dimensions, the w velocity must be supplied!");
      72          66 : }
      73             : 
      74             : ADReal
      75       94300 : WCNSFV2PInterfaceAreaSourceSink::computeQpResidual()
      76             : {
      77             : 
      78             :   // Useful Arguments
      79       94300 :   const auto state = determineState();
      80       94300 :   const auto elem_arg = makeElemArg(_current_elem);
      81       94300 :   const bool is_transient = _subproblem.isTransient();
      82             : 
      83             :   // Current Variables
      84       94300 :   const auto u = _u_var(elem_arg, state);
      85       94300 :   const auto rho_d = _rho_d(elem_arg, state);
      86       94300 :   const auto rho_d_grad = _rho_d.gradient(elem_arg, state);
      87       94300 :   const auto xi = _var(elem_arg, state);
      88       94300 :   const auto rho_m = _rho_mixture(elem_arg, state);
      89       94300 :   const auto f_d = _f_d(elem_arg, state);
      90       94300 :   const auto sigma = _sigma(elem_arg, state);
      91       94300 :   const auto rho_l = (rho_m - f_d * rho_d) / (1.0 - f_d + libMesh::TOLERANCE);
      92      188600 :   const auto complement_fd = std::max(_f_d_max - f_d, libMesh::TOLERANCE);
      93       94300 :   const auto f_d_o_xi = f_d / (_var(elem_arg, state) + libMesh::TOLERANCE) + libMesh::TOLERANCE;
      94             :   const auto f_d_o_xi_old =
      95      188600 :       f_d / (raw_value(_var(elem_arg, state)) + libMesh::TOLERANCE) + libMesh::TOLERANCE;
      96             : 
      97             :   // Adding bubble compressibility term
      98             :   ADReal material_time_derivative_rho_d = u * rho_d_grad(0);
      99       94300 :   if (_dim > 1)
     100             :   {
     101      188600 :     material_time_derivative_rho_d += (*_v_var)(elem_arg, state) * rho_d_grad(1);
     102       94300 :     if (_dim > 2)
     103             :     {
     104           0 :       material_time_derivative_rho_d += (*_w_var)(elem_arg, state) * rho_d_grad(2);
     105             :     }
     106             :   }
     107       94300 :   if (is_transient)
     108             :     material_time_derivative_rho_d +=
     109      134600 :         raw_value(_rho_d(elem_arg, state) - _rho_d(elem_arg, Moose::oldState())) / _dt;
     110      188600 :   const auto bubble_compressibility = material_time_derivative_rho_d * xi / 3.0;
     111             : 
     112             :   // Adding area growth due to added mass
     113             :   ADReal bubble_added_mass;
     114       94300 :   if (f_d < _cutoff_fraction)
     115           0 :     bubble_added_mass = raw_value(_rho_d(elem_arg, state)) *
     116           0 :                         (f_d * 6.0 / _particle_diameter(elem_arg, state) - _var(elem_arg, state));
     117             :   else
     118       94300 :     bubble_added_mass = 2. / 3. * _mass_exchange_coefficient(elem_arg, state) *
     119      188600 :                         (1.0 / (f_d + libMesh::TOLERANCE) - 2.0);
     120             : 
     121             :   // Model parameters
     122       94300 :   const auto db = _shape_factor * f_d_o_xi_old + libMesh::TOLERANCE;
     123             : 
     124      188600 :   ADRealVectorValue velocity(u);
     125       94300 :   if (_v_var)
     126       94300 :     velocity(1) = (*_v_var)(elem_arg, state);
     127       94300 :   if (_w_var)
     128           0 :     velocity(2) = (*_w_var)(elem_arg, state);
     129             : 
     130       94300 :   const ADReal velocity_norm = NS::computeSpeed<ADReal>(velocity);
     131      188600 :   const auto pressure_gradient = raw_value(_pressure.gradient(elem_arg, state));
     132             :   const Real pressure_grad_norm =
     133       94300 :       MooseUtils::isZero(pressure_gradient) ? 1e-42 : pressure_gradient.norm();
     134             : 
     135             :   const auto u_eps =
     136      188600 :       std::pow(velocity_norm * _characheristic_length(elem_arg, state) * pressure_grad_norm / rho_m,
     137      188600 :                1. / 3.);
     138             : 
     139             :   const auto interaction_prefactor =
     140      188600 :       Utility::pow<2>(f_d_o_xi) * u_eps / (std::pow(db, 11. / 3.) / complement_fd);
     141             : 
     142             :   // Adding coalescence term
     143       94300 :   const auto f_c = interaction_prefactor * _gamma_c * Utility::pow<2>(f_d);
     144      188600 :   const auto exp_c = std::exp(-_Kc * std::pow(db, 5. / 6.) * std::sqrt(rho_l / sigma) * u_eps);
     145             :   const auto s_rc = f_c * exp_c;
     146             : 
     147             :   // Adding breakage term
     148      188600 :   const auto f_b = interaction_prefactor * _gamma_b * f_d * (1. - f_d);
     149             :   const auto exp_b =
     150      282900 :       std::exp(-_Kb * sigma / (rho_l * std::pow(db, 5. / 3.) * Utility::pow<2>(u_eps)));
     151             :   const auto s_rb = f_b * exp_b;
     152             : 
     153      188600 :   return -bubble_added_mass + bubble_compressibility + s_rc - s_rb;
     154             : }