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 "PorousFlowFluidStateFlash.h"
      11             : 
      12             : InputParameters
      13        1406 : PorousFlowFluidStateFlash::validParams()
      14             : {
      15        1406 :   InputParameters params = PorousFlowFluidStateBase::validParams();
      16        1406 :   params.addClassDescription("Compositional flash calculations for use in fluid state classes");
      17        1406 :   return params;
      18           0 : }
      19             : 
      20         703 : PorousFlowFluidStateFlash::PorousFlowFluidStateFlash(const InputParameters & parameters)
      21         703 :   : PorousFlowFluidStateBase(parameters), _nr_max_its(42), _nr_tol(1.0e-12)
      22             : {
      23         703 : }
      24             : 
      25             : Real
      26          15 : PorousFlowFluidStateFlash::rachfordRice(Real x,
      27             :                                         std::vector<Real> & Zi,
      28             :                                         std::vector<Real> & Ki) const
      29             : {
      30             :   const std::size_t num_z = Zi.size();
      31             :   // Check that the sizes of the mass fractions and equilibrium constant vectors are correct
      32          15 :   if (Ki.size() != num_z + 1)
      33           0 :     mooseError("The number of mass fractions or equilibrium components passed to rachfordRice is "
      34             :                "not correct");
      35             : 
      36             :   Real f = 0.0;
      37             :   Real Z_total = 0.0;
      38             : 
      39          54 :   for (std::size_t i = 0; i < num_z; ++i)
      40             :   {
      41          39 :     f += Zi[i] * (Ki[i] - 1.0) / (1.0 + x * (Ki[i] - 1.0));
      42          39 :     Z_total += Zi[i];
      43             :   }
      44             : 
      45             :   // Add the last component (with total mass fraction = 1 - z_total)
      46          15 :   f += (1.0 - Z_total) * (Ki[num_z] - 1.0) / (1.0 + x * (Ki[num_z] - 1.0));
      47             : 
      48          15 :   return f;
      49             : }
      50             : 
      51             : Real
      52           6 : PorousFlowFluidStateFlash::rachfordRiceDeriv(Real x,
      53             :                                              std::vector<Real> & Zi,
      54             :                                              std::vector<Real> & Ki) const
      55             : {
      56             :   const std::size_t num_Z = Zi.size();
      57             :   // Check that the sizes of the mass fractions and equilibrium constant vectors are correct
      58           6 :   if (Ki.size() != num_Z + 1)
      59           0 :     mooseError("The number of mass fractions or equilibrium components passed to rachfordRice is "
      60             :                "not correct");
      61             : 
      62             :   Real df = 0.0;
      63             :   Real Z_total = 0.0;
      64             : 
      65          22 :   for (std::size_t i = 0; i < num_Z; ++i)
      66             :   {
      67          16 :     df -= Zi[i] * (Ki[i] - 1.0) * (Ki[i] - 1.0) / (1.0 + x * (Ki[i] - 1.0)) /
      68             :           (1.0 + x * (Ki[i] - 1.0));
      69          16 :     Z_total += Zi[i];
      70             :   }
      71             : 
      72             :   // Add the last component (with total mass fraction = 1 - z_total)
      73           6 :   df -= (1.0 - Z_total) * (Ki[num_Z] - 1.0) * (Ki[num_Z] - 1.0) / (1.0 + x * (Ki[num_Z] - 1.0)) /
      74             :         (1.0 + x * (Ki[num_Z] - 1.0));
      75             : 
      76           6 :   return df;
      77             : }
      78             : 
      79             : Real
      80           1 : PorousFlowFluidStateFlash::vaporMassFraction(Real Z0, Real K0, Real K1) const
      81             : {
      82           1 :   return (Z0 * (K1 - K0) - (K1 - 1.0)) / ((K0 - 1.0) * (K1 - 1.0));
      83             : }
      84             : 
      85             : ADReal
      86       86764 : PorousFlowFluidStateFlash::vaporMassFraction(const ADReal & Z0,
      87             :                                              const ADReal & K0,
      88             :                                              const ADReal & K1) const
      89             : {
      90       86764 :   return (Z0 * (K1 - K0) - (K1 - 1.0)) / ((K0 - 1.0) * (K1 - 1.0));
      91             : }
      92             : 
      93             : Real
      94           2 : PorousFlowFluidStateFlash::vaporMassFraction(std::vector<Real> & Zi, std::vector<Real> & Ki) const
      95             : {
      96             :   // Check that the sizes of the mass fractions and equilibrium constant vectors are correct
      97           2 :   if (Ki.size() != Zi.size() + 1)
      98           0 :     mooseError("The number of mass fractions or equilibrium components passed to rachfordRice is "
      99             :                "not correct");
     100             :   Real v;
     101             : 
     102             :   // If there are only two components, an analytical solution is possible
     103           2 :   if (Ki.size() == 2)
     104           1 :     v = vaporMassFraction(Zi[0], Ki[0], Ki[1]);
     105             :   else
     106             :   {
     107             :     // More than two components - solve the Rachford-Rice equation using
     108             :     // Newton-Raphson method
     109             :     // Initial guess for vapor mass fraction
     110             :     Real v0 = 0.5;
     111             :     unsigned int iter = 0;
     112             : 
     113           5 :     while (std::abs(rachfordRice(v0, Zi, Ki)) > _nr_tol)
     114             :     {
     115           4 :       v0 = v0 - rachfordRice(v0, Zi, Ki) / rachfordRiceDeriv(v0, Zi, Ki);
     116           4 :       iter++;
     117             : 
     118           4 :       if (iter > _nr_max_its)
     119             :         break;
     120             :     }
     121             :     v = v0;
     122             :   }
     123           2 :   return v;
     124             : }