GCC Code Coverage Report

Directory:	../../../builds/dumux-repositories/
File:	/builds/dumux-repositories/dumux/examples/biomineralization/material/fluidsystems/icpcomplexsalinitybrine.hh
Date:	2024-09-21 20:52:54

	Exec	Total	Coverage
Lines:	37	38	97.4%
Functions:	3	3	100.0%
Branches:	8	32	25.0%

  
      Line
      Branch
      Exec
      Source
    
      // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
    
      // vi: set et ts=4 sw=4 sts=4:
    
      //
    
      // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
    
      // SPDX-License-Identifier: GPL-3.0-or-later
    
      //
    
      #ifndef DUMUX_ICP_COMPLEX_SALINITY_BRINE_FLUID_SYSTEM_HH
    
      #define DUMUX_ICP_COMPLEX_SALINITY_BRINE_FLUID_SYSTEM_HH
    
      // ## The brine adapter (`icpcomplexsalinitybrine.hh`)
    
      //
    
      // This file contains the brineadapter class__ which defines
    
      // how to adapt values for communication between
    
      // the "full" fluidsystem accounting for 4 components influencing the brine properties
    
      // and the "brine" fluidsystem assuming water and a single salt determining the brine properties.
    
      //
    
      // [[content]]
    
      //
    
      // ### Include files
    
      // [[details]] includes
    
      // we include all necessary components
    
      #include <dune/common/math.hh>
    
      #include <dumux/material/fluidsystems/base.hh>
    
      #include <dumux/material/constants.hh>
    
      #include <dumux/material/components/h2o.hh>
    
      #include <dumux/material/components/sodiumion.hh>
    
      #include <dumux/material/components/chlorideion.hh>
    
      #include <dumux/material/components/calciumion.hh>
    
      #include <dumux/material/components/tabulatedcomponent.hh>
    
      #include <dumux/common/exceptions.hh>
    
      #include <dumux/io/name.hh>
    
      // [[/details]]
    
      //
    
      // ### The brine adapter class
    
      // In ICPComplexSalinityBrine we make the transition from 3 components additional to water contributing to salinity to the single component that is assumed in the brine fluid system.
    
      // We enter the namespace Dumux. All Dumux functions and classes are in a namespace Dumux, to make sure they don`t clash with symbols from other libraries you may want to use in conjunction with Dumux.
    
      // [[codeblock]]
    
      namespace Dumux::FluidSystems {
    
      template< class Scalar, class H2OType = Components::TabulatedComponent<Dumux::Components::H2O<Scalar>> >
    
      class ICPComplexSalinityBrine : public Base< Scalar, ICPComplexSalinityBrine<Scalar, H2OType>>
    
      {
    
          using ThisType = ICPComplexSalinityBrine<Scalar, H2OType>;
    
          using Base = Dumux::FluidSystems::Base<Scalar, ThisType>;
    
      public:
    
          //! export the involved components
    
          using H2O = H2OType;
    
          using Na = Components::SodiumIon<Scalar>;
    
          using Cl = Components::ChlorideIon<Scalar>;
    
          using Ca = Components::CalciumIon<Scalar>;
    
          // [[/codeblock]]
    
          //
    
          // #### Component and phase definitions
    
          // With the following function we define what phases and components will be used by the fluid system and define the indices used to distinguish  phases and components in the course of the simulation.
    
          // [[codeblock]]
    
          // We use convenient declarations that we derive from the property system.
    
          static constexpr int numPhases = 1;     //!< Number of phases in the fluid system
    
          static constexpr int numComponents = 4; //!< Number of components in the fluid system (H2O, Na, Cl, Ca)
    
          static constexpr int phase0Idx = 0; //!< Index of the first (and only) phase
    
          static constexpr int liquidPhaseIdx = phase0Idx; //!< The one considered phase is liquid
    
          static constexpr int H2OIdx = 0;         //!< index of the water component
    
          static constexpr int NaIdx = 1;        //!< index of the Na component
    
          static constexpr int ClIdx = 2;        //!< index of the Cl component
    
          static constexpr int CaIdx = 3;        //!< index of the Ca component
    
          static constexpr int comp0Idx = H2OIdx;  //!< index of the first component
    
          static constexpr int comp1Idx = NaIdx; //!< index of the second component
    
          static constexpr int comp2Idx = ClIdx; //!< index of the third component
    
          static constexpr int comp3Idx = CaIdx; //!< index of the fourth component
    
          // the fluid phase index
    
          static std::string phaseName(int phaseIdx = liquidPhaseIdx)
    
          {
    
              assert(phaseIdx == liquidPhaseIdx);
    
              return IOName::liquidPhase();
    
          }
    
          // the miscibility
    
          static constexpr bool isMiscible()
    
          {
    
              return false;
    
          }
    
          // we do not have a gas phase
    
          static constexpr bool isGas(int phaseIdx = liquidPhaseIdx)
    
          {
    
              assert(phaseIdx == liquidPhaseIdx);
    
              return false;
    
          }
    
          // We need a ideal mixture
    
          static constexpr bool isIdealMixture(int phaseIdx = liquidPhaseIdx)
    
          {
    
              assert(phaseIdx == liquidPhaseIdx);
    
              return true;
    
          }
    
          // phase compressibility
    
          static constexpr bool isCompressible(int phaseIdx = liquidPhaseIdx)
    
          {
    
              assert(phaseIdx == liquidPhaseIdx);
    
              return H2O::liquidIsCompressible();
    
          }
    
          // no gas, thus no ideal gas
    
          static constexpr bool isIdealGas(int phaseIdx = liquidPhaseIdx)
    
          {
    
              assert(phaseIdx == liquidPhaseIdx);
    
              return false; /*we're a liquid!*/
    
          }
    
      // [[/codeblock]]
    
      // [[exclude]]
    
          // The component names
    
          static std::string componentName(int compIdx)
    
          {
    
              switch (compIdx)
    
              {
    
                  case H2OIdx: return H2O::name();
    
                  case CaIdx: return Ca::name();
    
                  case NaIdx: return Na::name();
    
                  case ClIdx: return Cl::name();
    
              };
    
              DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
    
          }
    
          // The component molar masses
    
        1/2✓ Branch 0 taken 20766176 times.
✗ Branch 1 not taken.

      20766176
          static Scalar molarMass(int compIdx)
    
          {
    
              switch (compIdx)
    
              {
    
                  case H2OIdx: return H2O::molarMass();
    
                  case CaIdx: return Ca::molarMass();
    
                  case NaIdx: return Na::molarMass();
    
                  case ClIdx: return Cl::molarMass();
    
              };
    
      ✗
              DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
    
          }
    
          //initializing
    
           static void init()
    
           {
    
               init(/*tempMin=*/273.15,
    
                    /*tempMax=*/623.15,
    
                    /*numTemp=*/100,
    
                    /*pMin=*/-10.,
    
                    /*pMax=*/20e6,
    
                    /*numP=*/200);
    
           }
    
          static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
    
                           Scalar pressMin, Scalar pressMax, unsigned nPress)
    
          {
    
              if (H2O::isTabulated)
    
              {
    
                  std::cout << "Initializing tables for the H2O fluid properties ("
    
                            << nTemp*nPress
    
                            << " entries).\n";
    
                  H2O::init(tempMin, tempMax, nTemp, pressMin, pressMax, nPress);
    
              }
    
          }
    
      // [[/exclude]]
    
          //
    
          // #### The Component-Phase Interactions
    
          // In the following, we define the component-phase interactions such as // each component's fugacity coefficient in brine
    
          // and each component's binary diffusion coefficient with water in brine
    
          // [[codeblock]]
    
          // The component fugacity coefficients
    
          template <class FluidState>
    
          static Scalar fugacityCoefficient(const FluidState &fluidState,
    
                                            int phaseIdx,
    
                                            int compIdx)
    
          {
    
              assert(0 <= phaseIdx  && phaseIdx < numPhases);
    
              assert(0 <= compIdx  && compIdx < numComponents);
    
              if (phaseIdx == compIdx)
    
                  // We could calculate the real fugacity coefficient of
    
                  // the component in the fluid. Probably that's not worth
    
                  // the effort, since the fugacity coefficient of the other
    
                  // component is infinite anyway...
    
                  return 1.0;
    
              return std::numeric_limits<Scalar>::infinity();
    
          }
    
          // The binary diffusion coefficients of the components and the phases main component
    
          template <class FluidState>
    
          static Scalar binaryDiffusionCoefficient(const FluidState& fluidState,
    
                                                   int phaseIdx,
    
                                                   int compIIdx,
    
                                                   int compJIdx)
    
          {
    
              if (phaseIdx == liquidPhaseIdx)
    
              {
    
                  if (compIIdx > compJIdx)
    
                  {
    
                      using std::swap;
    
                      swap(compIIdx, compJIdx);
    
                  }
    
                  if (compJIdx == NaIdx || compJIdx == ClIdx || compJIdx == CaIdx)
    
                      return 1.587e-9;  //[m²/s]    //J. Phys. D: Appl. Phys. 40 (2007) 2769-2776
    
                  else
    
                      DUNE_THROW(Dune::NotImplemented, "Binary diffusion coefficient of components "
    
                                                       << compIIdx << " and " << compJIdx
    
                                                       << " in phase " << phaseIdx);
    
               }
    
               DUNE_THROW(Dune::InvalidStateException, "Invalid phase index: " << phaseIdx);
    
          }
    
          // [[/codeblock]]
    
          //
    
          // #### The Fluid (Brine) Properties
    
          // In the following, all functions defining the phase properties are given:
    
          // the density, viscosity, enthalpy, thermal conductivities, and heat capacities
    
          // of each phase depending on temperature, pressure, and phase composition
    
          // [[codeblock]]
    
          // The phase density
    
          template <class FluidState>
    
      5191544
          static Scalar density(const FluidState& fluidState, int phaseIdx = liquidPhaseIdx)
    
          {
    
        1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5191544 times.

      5191544
              assert(phaseIdx == liquidPhaseIdx);
    
      5191544
              const Scalar temperature = fluidState.temperature(phaseIdx);
    
      5191544
              const Scalar pressure = fluidState.pressure(phaseIdx);
    
      5191544
              const Scalar xNaCl = fluidState.massFraction(phaseIdx, NaIdx)
    
      5191544
                                 + fluidState.massFraction(phaseIdx, ClIdx)
    
      5191544
                                 + fluidState.massFraction(phaseIdx, CaIdx);
    
              using std::max;
    
      5191544
              const Scalar TempC = temperature - 273.15;
    
      5191544
              const Scalar pMPa = pressure/1.0E6;
    
        1/2✓ Branch 0 taken 5191544 times.
✗ Branch 1 not taken.

      5191544
              const Scalar salinity = max(0.0, xNaCl);
    
      5191544
              const Scalar rhow = H2O::liquidDensity(temperature, pressure);
    
      15574632
              const Scalar density = rhow + 1000*salinity*(0.668 +
    
      10383088
                                                           0.44*salinity +
    
      15574632
                                                           1.0E-6*(300*pMPa -
    
      10383088
                                                                   2400*pMPa*salinity +
    
      10383088
                                                                   TempC*(80.0 +
    
      10383088
                                                                          3*TempC -
    
      10383088
                                                                          3300*salinity -
    
      10383088
                                                                          13*pMPa +
    
      5191544
                                                                          47*pMPa*salinity)));
    
        1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5191544 times.

      5191544
              assert(density > 0.0);
    
      5191544
              return density;
    
          }
    
          // The phase viscosity
    
          template <class FluidState>
    
      2595772
          static Scalar viscosity(const FluidState& fluidState, int phaseIdx = liquidPhaseIdx)
    
          {
    
        1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2595772 times.

      2595772
              assert(phaseIdx == liquidPhaseIdx);
    
      2595772
              const Scalar temperature = fluidState.temperature(phaseIdx);
    
      2595772
              const Scalar xNaCl = fluidState.massFraction(phaseIdx, NaIdx)
    
      2595772
                                 + fluidState.massFraction(phaseIdx, ClIdx)
    
      2595772
                                 + fluidState.massFraction(phaseIdx, CaIdx);
    
              using std::pow;
    
              using Dune::power;
    
              using std::exp;
    
              using std::max;
    
        1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2595772 times.

      2595772
              const Scalar T = max(temperature, 275.0);
    
        1/2✓ Branch 0 taken 2595772 times.
✗ Branch 1 not taken.

      2595772
              const Scalar salinity = max(0.0, xNaCl);
    
      2595772
              const Scalar T_C = T - 273.15;
    
      2595772
              const Scalar A = ((0.42*power((pow(salinity, 0.8)-0.17), 2)) + 0.045)*pow(T_C, 0.8);
    
      2595772
              const Scalar mu_brine = 0.1 + (0.333*salinity) + (1.65+(91.9*salinity*salinity*salinity))*exp(-A); // [cP]
    
        1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2595772 times.

      2595772
              assert(mu_brine > 0.0);
    
      2595772
              return mu_brine/1000.0; // [Pa·s]
    
          }
    
      // [[/codeblock]]
    
      // [[exclude]]
    
          // The vapor pressure
    
          template <class FluidState>
    
          static Scalar vaporPressure(const FluidState& fluidState, int compIdx)
    
          {
    
              if (compIdx == H2OIdx)
    
              {
    
                  // simplified version of Eq 2.29 in Vishal Jambhekar's Promo
    
                  const Scalar temperature = fluidState.temperature(H2OIdx);
    
                  return H2O::vaporPressure(temperature)/fluidState.massFraction(phase0Idx, H2OIdx);
    
              }
    
              else if (compIdx == NaIdx)
    
                  DUNE_THROW(Dune::NotImplemented, "Na::vaporPressure(t)");
    
              else if (compIdx == ClIdx)
    
                  DUNE_THROW(Dune::NotImplemented, "Cl::vaporPressure(t)");
    
              else if (compIdx == CaIdx)
    
                  DUNE_THROW(Dune::NotImplemented, "Ca::vaporPressure(t)");
    
              else
    
                  DUNE_THROW(Dune::NotImplemented, "Invalid component index " << compIdx);
    
          }
    
          // The phase enthalpies
    
          template <class FluidState>
    
          static Scalar enthalpy(const FluidState& fluidState, int phaseIdx)
    
          {
    
              /*Numerical coefficients from PALLISER*/
    
              static constexpr Scalar f[] = { 2.63500E-1, 7.48368E-6, 1.44611E-6, -3.80860E-10 };
    
              /*Numerical coefficients from MICHAELIDES for the enthalpy of brine*/
    
              static constexpr Scalar a[4][3] = { { +9633.6, -4080.0, +286.49 },
    
                                                  { +166.58, +68.577, -4.6856 },
    
                                                  { -0.90963, -0.36524, +0.249667E-1 },
    
                                                  { +0.17965E-2, +0.71924E-3, -0.4900E-4 } };
    
              const Scalar p = fluidState.pressure(phaseIdx);
    
              const Scalar T = fluidState.temperature(phaseIdx);
    
              const Scalar theta = T - 273.15;
    
              const Scalar salSat = f[0] + f[1]*theta + f[2]*theta*theta + f[3]*theta*theta*theta;
    
              /*Regularization*/
    
              using std::min;
    
              using std::max;
    
              const Scalar xNaCl = fluidState.massFraction(phaseIdx, NaIdx)
    
                                 + fluidState.massFraction(phaseIdx, ClIdx)
    
                                 + fluidState.massFraction(phaseIdx, CaIdx);
    
              const Scalar salinity = min(max(xNaCl,0.0), salSat);
    
              const Scalar hw = H2O::liquidEnthalpy(T, p)/1E3; /* kJ/kg */
    
              /*DAUBERT and DANNER*/
    
              const Scalar h_NaCl = (3.6710E4*T + 0.5*(6.2770E1)*T*T - ((6.6670E-2)/3)*T*T*T
    
                                    + ((2.8000E-5)/4)*(T*T*T*T))/(58.44E3)- 2.045698e+02; /* U [kJ/kg] */
    
              const Scalar m = (1E3/58.44)*(salinity/(1-salinity));
    
              using Dune::power;
    
              Scalar d_h = 0;
    
              for (int i = 0; i<=3; i++)
    
                  for (int j=0; j<=2; j++)
    
                      d_h = d_h + a[i][j] * power(theta, i) * power(m, j);
    
              /* heat of dissolution for halite according to Michaelides 1971 */
    
              const Scalar delta_h = (4.184/(1E3 + (58.44 * m)))*d_h;
    
              /* Enthalpy of brine without any dissolved gas */
    
              const Scalar h_ls1 =(1-salinity)*hw + salinity*h_NaCl + salinity*delta_h; /* kJ/kg */
    
              return h_ls1*1E3; /*J/kg*/
    
          }
    
          // The molar density
    
          template <class FluidState>
    
          static Scalar molarDensity(const FluidState& fluidState, int phaseIdx = liquidPhaseIdx)
    
          {
    
      5191544
              return density(fluidState, phaseIdx)/fluidState.averageMolarMass(phaseIdx);
    
          }
    
          template <class FluidState>
    
          static Scalar diffusionCoefficient(const FluidState& fluidState, int phaseIdx, int compIdx)
    
          {
    
              DUNE_THROW(Dune::NotImplemented, "FluidSystems::ICPComplexSalinityBrine::diffusionCoefficient()");
    
          }
    
          // The thermal conductivity
    
          template <class FluidState>
    
          static Scalar thermalConductivity(const FluidState& fluidState, int phaseIdx)
    
          {
    
              if (phaseIdx == liquidPhaseIdx)
    
                  return H2O::liquidThermalConductivity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
    
              DUNE_THROW(Dune::InvalidStateException, "Invalid phase index: " << phaseIdx);
    
          }
    
          // The phase heat capacity
    
          template <class FluidState>
    
          static Scalar heatCapacity(const FluidState &fluidState, int phaseIdx)
    
          {
    
              if (phaseIdx == liquidPhaseIdx)
    
                  return H2O::liquidHeatCapacity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
    
              DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    
          }
    
      };
    
      } // end namespace Dumux::FluidSystems
    
      // [[/exclude]]
    
      // [[/content]]
    
      #endif

Line	Branch	Exec	Source
1			// -- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 --
2			// vi: set et ts=4 sw=4 sts=4:
3			//
4			// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5			// SPDX-License-Identifier: GPL-3.0-or-later
6			//
7
8			#ifndef DUMUX_ICP_COMPLEX_SALINITY_BRINE_FLUID_SYSTEM_HH
9			#define DUMUX_ICP_COMPLEX_SALINITY_BRINE_FLUID_SYSTEM_HH
10
11			// ## The brine adapter (`icpcomplexsalinitybrine.hh`)
12			//
13			// This file contains the brineadapter class__ which defines
14			// how to adapt values for communication between
15			// the "full" fluidsystem accounting for 4 components influencing the brine properties
16			// and the "brine" fluidsystem assuming water and a single salt determining the brine properties.
17			//
18			// [[content]]
19			//
20			// ### Include files
21			// [[details]] includes
22			// we include all necessary components
23			#include <dune/common/math.hh>
24
25			#include <dumux/material/fluidsystems/base.hh>
26			#include <dumux/material/constants.hh>
27			#include <dumux/material/components/h2o.hh>
28			#include <dumux/material/components/sodiumion.hh>
29			#include <dumux/material/components/chlorideion.hh>
30			#include <dumux/material/components/calciumion.hh>
31			#include <dumux/material/components/tabulatedcomponent.hh>
32
33			#include <dumux/common/exceptions.hh>
34
35			#include <dumux/io/name.hh>
36			// [[/details]]
37			//
38			// ### The brine adapter class
39			// In ICPComplexSalinityBrine we make the transition from 3 components additional to water contributing to salinity to the single component that is assumed in the brine fluid system.
40			// We enter the namespace Dumux. All Dumux functions and classes are in a namespace Dumux, to make sure they don`t clash with symbols from other libraries you may want to use in conjunction with Dumux.
41			// [[codeblock]]
42			namespace Dumux::FluidSystems {
43
44			template< class Scalar, class H2OType = Components::TabulatedComponent<Dumux::Components::H2O<Scalar>> >
45			class ICPComplexSalinityBrine : public Base< Scalar, ICPComplexSalinityBrine<Scalar, H2OType>>
46			{
47			using ThisType = ICPComplexSalinityBrine<Scalar, H2OType>;
48			using Base = Dumux::FluidSystems::Base<Scalar, ThisType>;
49
50			public:
51			//! export the involved components
52			using H2O = H2OType;
53			using Na = Components::SodiumIon<Scalar>;
54			using Cl = Components::ChlorideIon<Scalar>;
55			using Ca = Components::CalciumIon<Scalar>;
56			// [[/codeblock]]
57			//
58			// #### Component and phase definitions
59			// With the following function we define what phases and components will be used by the fluid system and define the indices used to distinguish phases and components in the course of the simulation.
60			// [[codeblock]]
61			// We use convenient declarations that we derive from the property system.
62			static constexpr int numPhases = 1; //!< Number of phases in the fluid system
63			static constexpr int numComponents = 4; //!< Number of components in the fluid system (H2O, Na, Cl, Ca)
64
65			static constexpr int phase0Idx = 0; //!< Index of the first (and only) phase
66			static constexpr int liquidPhaseIdx = phase0Idx; //!< The one considered phase is liquid
67
68			static constexpr int H2OIdx = 0; //!< index of the water component
69			static constexpr int NaIdx = 1; //!< index of the Na component
70			static constexpr int ClIdx = 2; //!< index of the Cl component
71			static constexpr int CaIdx = 3; //!< index of the Ca component
72			static constexpr int comp0Idx = H2OIdx; //!< index of the first component
73			static constexpr int comp1Idx = NaIdx; //!< index of the second component
74			static constexpr int comp2Idx = ClIdx; //!< index of the third component
75			static constexpr int comp3Idx = CaIdx; //!< index of the fourth component
76
77			// the fluid phase index
78			static std::string phaseName(int phaseIdx = liquidPhaseIdx)
79			{
80			assert(phaseIdx == liquidPhaseIdx);
81			return IOName::liquidPhase();
82			}
83
84			// the miscibility
85			static constexpr bool isMiscible()
86			{
87			return false;
88			}
89
90			// we do not have a gas phase
91			static constexpr bool isGas(int phaseIdx = liquidPhaseIdx)
92			{
93			assert(phaseIdx == liquidPhaseIdx);
94			return false;
95			}
96
97			// We need a ideal mixture
98			static constexpr bool isIdealMixture(int phaseIdx = liquidPhaseIdx)
99			{
100			assert(phaseIdx == liquidPhaseIdx);
101			return true;
102			}
103
104			// phase compressibility
105			static constexpr bool isCompressible(int phaseIdx = liquidPhaseIdx)
106			{
107			assert(phaseIdx == liquidPhaseIdx);
108			return H2O::liquidIsCompressible();
109			}
110
111			// no gas, thus no ideal gas
112			static constexpr bool isIdealGas(int phaseIdx = liquidPhaseIdx)
113			{
114			assert(phaseIdx == liquidPhaseIdx);
115			return false; /we're a liquid!/
116			}
117			// [[/codeblock]]
118			// [[exclude]]
119			// The component names
120			static std::string componentName(int compIdx)
121			{
122			switch (compIdx)
123			{
124			case H2OIdx: return H2O::name();
125			case CaIdx: return Ca::name();
126			case NaIdx: return Na::name();
127			case ClIdx: return Cl::name();
128			};
129			DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
130			}
131
132			// The component molar masses
133	1/2 ✓ Branch 0 taken 20766176 times. ✗ Branch 1 not taken.	20766176	static Scalar molarMass(int compIdx)
134			{
135			switch (compIdx)
136			{
137			case H2OIdx: return H2O::molarMass();
138			case CaIdx: return Ca::molarMass();
139			case NaIdx: return Na::molarMass();
140			case ClIdx: return Cl::molarMass();
141			};
142		✗	DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
143			}
144
145			//initializing
146			static void init()
147			{
148			init(/tempMin=/273.15,
149			/tempMax=/623.15,
150			/numTemp=/100,
151			/pMin=/-10.,
152			/pMax=/20e6,
153			/numP=/200);
154			}
155			static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
156			Scalar pressMin, Scalar pressMax, unsigned nPress)
157			{
158
159			if (H2O::isTabulated)
160			{
161			std::cout << "Initializing tables for the H2O fluid properties ("
162			<< nTemp*nPress
163			<< " entries).\n";
164
165			H2O::init(tempMin, tempMax, nTemp, pressMin, pressMax, nPress);
166			}
167			}
168			// [[/exclude]]
169			//
170			// #### The Component-Phase Interactions
171			// In the following, we define the component-phase interactions such as // each component's fugacity coefficient in brine
172			// and each component's binary diffusion coefficient with water in brine
173			// [[codeblock]]
174			// The component fugacity coefficients
175			template <class FluidState>
176			static Scalar fugacityCoefficient(const FluidState &fluidState,
177			int phaseIdx,
178			int compIdx)
179			{
180			assert(0 <= phaseIdx && phaseIdx < numPhases);
181			assert(0 <= compIdx && compIdx < numComponents);
182
183			if (phaseIdx == compIdx)
184			// We could calculate the real fugacity coefficient of
185			// the component in the fluid. Probably that's not worth
186			// the effort, since the fugacity coefficient of the other
187			// component is infinite anyway...
188			return 1.0;
189			return std::numeric_limits<Scalar>::infinity();
190			}
191
192			// The binary diffusion coefficients of the components and the phases main component
193			template <class FluidState>
194			static Scalar binaryDiffusionCoefficient(const FluidState& fluidState,
195			int phaseIdx,
196			int compIIdx,
197			int compJIdx)
198			{
199			if (phaseIdx == liquidPhaseIdx)
200			{
201			if (compIIdx > compJIdx)
202			{
203			using std::swap;
204			swap(compIIdx, compJIdx);
205			}
206			if (compJIdx == NaIdx \|\| compJIdx == ClIdx \|\| compJIdx == CaIdx)
207			return 1.587e-9; //[m²/s] //J. Phys. D: Appl. Phys. 40 (2007) 2769-2776
208			else
209			DUNE_THROW(Dune::NotImplemented, "Binary diffusion coefficient of components "
210			<< compIIdx << " and " << compJIdx
211			<< " in phase " << phaseIdx);
212			}
213
214			DUNE_THROW(Dune::InvalidStateException, "Invalid phase index: " << phaseIdx);
215			}
216			// [[/codeblock]]
217			//
218			// #### The Fluid (Brine) Properties
219			// In the following, all functions defining the phase properties are given:
220			// the density, viscosity, enthalpy, thermal conductivities, and heat capacities
221			// of each phase depending on temperature, pressure, and phase composition
222			// [[codeblock]]
223			// The phase density
224			template <class FluidState>
225		5191544	static Scalar density(const FluidState& fluidState, int phaseIdx = liquidPhaseIdx)
226			{
227	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 5191544 times.	5191544	assert(phaseIdx == liquidPhaseIdx);
228		5191544	const Scalar temperature = fluidState.temperature(phaseIdx);
229		5191544	const Scalar pressure = fluidState.pressure(phaseIdx);
230		5191544	const Scalar xNaCl = fluidState.massFraction(phaseIdx, NaIdx)
231		5191544	+ fluidState.massFraction(phaseIdx, ClIdx)
232		5191544	+ fluidState.massFraction(phaseIdx, CaIdx);
233
234			using std::max;
235		5191544	const Scalar TempC = temperature - 273.15;
236		5191544	const Scalar pMPa = pressure/1.0E6;
237	1/2 ✓ Branch 0 taken 5191544 times. ✗ Branch 1 not taken.	5191544	const Scalar salinity = max(0.0, xNaCl);
238
239		5191544	const Scalar rhow = H2O::liquidDensity(temperature, pressure);
240		15574632	const Scalar density = rhow + 1000salinity(0.668 +
241		10383088	0.44*salinity +
242		15574632	1.0E-6(300pMPa -
243		10383088	2400pMPasalinity +
244		10383088	TempC*(80.0 +
245		10383088	3*TempC -
246		10383088	3300*salinity -
247		10383088	13*pMPa +
248		5191544	47pMPasalinity)));
249	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 5191544 times.	5191544	assert(density > 0.0);
250		5191544	return density;
251			}
252
253
254
255			// The phase viscosity
256			template <class FluidState>
257		2595772	static Scalar viscosity(const FluidState& fluidState, int phaseIdx = liquidPhaseIdx)
258			{
259	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 2595772 times.	2595772	assert(phaseIdx == liquidPhaseIdx);
260		2595772	const Scalar temperature = fluidState.temperature(phaseIdx);
261		2595772	const Scalar xNaCl = fluidState.massFraction(phaseIdx, NaIdx)
262		2595772	+ fluidState.massFraction(phaseIdx, ClIdx)
263		2595772	+ fluidState.massFraction(phaseIdx, CaIdx);
264
265			using std::pow;
266			using Dune::power;
267			using std::exp;
268			using std::max;
269	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 2595772 times.	2595772	const Scalar T = max(temperature, 275.0);
270	1/2 ✓ Branch 0 taken 2595772 times. ✗ Branch 1 not taken.	2595772	const Scalar salinity = max(0.0, xNaCl);
271
272		2595772	const Scalar T_C = T - 273.15;
273		2595772	const Scalar A = ((0.42power((pow(salinity, 0.8)-0.17), 2)) + 0.045)pow(T_C, 0.8);
274		2595772	const Scalar mu_brine = 0.1 + (0.333salinity) + (1.65+(91.9salinitysalinitysalinity))*exp(-A); // [cP]
275	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 2595772 times.	2595772	assert(mu_brine > 0.0);
276		2595772	return mu_brine/1000.0; // [Pa·s]
277			}
278			// [[/codeblock]]
279			// [[exclude]]
280			// The vapor pressure
281			template <class FluidState>
282			static Scalar vaporPressure(const FluidState& fluidState, int compIdx)
283			{
284			if (compIdx == H2OIdx)
285			{
286			// simplified version of Eq 2.29 in Vishal Jambhekar's Promo
287			const Scalar temperature = fluidState.temperature(H2OIdx);
288			return H2O::vaporPressure(temperature)/fluidState.massFraction(phase0Idx, H2OIdx);
289			}
290			else if (compIdx == NaIdx)
291			DUNE_THROW(Dune::NotImplemented, "Na::vaporPressure(t)");
292			else if (compIdx == ClIdx)
293			DUNE_THROW(Dune::NotImplemented, "Cl::vaporPressure(t)");
294			else if (compIdx == CaIdx)
295			DUNE_THROW(Dune::NotImplemented, "Ca::vaporPressure(t)");
296			else
297			DUNE_THROW(Dune::NotImplemented, "Invalid component index " << compIdx);
298			}
299
300			// The phase enthalpies
301			template <class FluidState>
302			static Scalar enthalpy(const FluidState& fluidState, int phaseIdx)
303			{
304			/Numerical coefficients from PALLISER/
305			static constexpr Scalar f[] = { 2.63500E-1, 7.48368E-6, 1.44611E-6, -3.80860E-10 };
306
307			/Numerical coefficients from MICHAELIDES for the enthalpy of brine/
308			static constexpr Scalar a[4][3] = { { +9633.6, -4080.0, +286.49 },
309			{ +166.58, +68.577, -4.6856 },
310			{ -0.90963, -0.36524, +0.249667E-1 },
311			{ +0.17965E-2, +0.71924E-3, -0.4900E-4 } };
312
313			const Scalar p = fluidState.pressure(phaseIdx);
314			const Scalar T = fluidState.temperature(phaseIdx);
315			const Scalar theta = T - 273.15;
316			const Scalar salSat = f[0] + f[1]theta + f[2]thetatheta + f[3]thetathetatheta;
317
318			/Regularization/
319			using std::min;
320			using std::max;
321			const Scalar xNaCl = fluidState.massFraction(phaseIdx, NaIdx)
322			+ fluidState.massFraction(phaseIdx, ClIdx)
323			+ fluidState.massFraction(phaseIdx, CaIdx);
324			const Scalar salinity = min(max(xNaCl,0.0), salSat);
325
326			const Scalar hw = H2O::liquidEnthalpy(T, p)/1E3; /* kJ/kg */
327
328			/DAUBERT and DANNER/
329			const Scalar h_NaCl = (3.6710E4T + 0.5(6.2770E1)TT - ((6.6670E-2)/3)TT*T
330			+ ((2.8000E-5)/4)(TTTT))/(58.44E3)- 2.045698e+02; /* U [kJ/kg] */
331
332			const Scalar m = (1E3/58.44)*(salinity/(1-salinity));
333
334			using Dune::power;
335			Scalar d_h = 0;
336			for (int i = 0; i<=3; i++)
337			for (int j=0; j<=2; j++)
338			d_h = d_h + a[i][j] * power(theta, i) * power(m, j);
339
340			/* heat of dissolution for halite according to Michaelides 1971 */
341			const Scalar delta_h = (4.184/(1E3 + (58.44 * m)))*d_h;
342
343			/* Enthalpy of brine without any dissolved gas */
344			const Scalar h_ls1 =(1-salinity)hw + salinityh_NaCl + salinitydelta_h; / kJ/kg */
345			return h_ls11E3; /J/kg*/
346			}
347			// The molar density
348			template <class FluidState>
349			static Scalar molarDensity(const FluidState& fluidState, int phaseIdx = liquidPhaseIdx)
350			{
351		5191544	return density(fluidState, phaseIdx)/fluidState.averageMolarMass(phaseIdx);
352			}
353
354			template <class FluidState>
355			static Scalar diffusionCoefficient(const FluidState& fluidState, int phaseIdx, int compIdx)
356			{
357			DUNE_THROW(Dune::NotImplemented, "FluidSystems::ICPComplexSalinityBrine::diffusionCoefficient()");
358			}
359
360			// The thermal conductivity
361			template <class FluidState>
362			static Scalar thermalConductivity(const FluidState& fluidState, int phaseIdx)
363			{
364			if (phaseIdx == liquidPhaseIdx)
365			return H2O::liquidThermalConductivity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
366			DUNE_THROW(Dune::InvalidStateException, "Invalid phase index: " << phaseIdx);
367			}
368
369			// The phase heat capacity
370			template <class FluidState>
371			static Scalar heatCapacity(const FluidState &fluidState, int phaseIdx)
372			{
373			if (phaseIdx == liquidPhaseIdx)
374			return H2O::liquidHeatCapacity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
375			DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
376			}
377			};
378
379			} // end namespace Dumux::FluidSystems
380			// [[/exclude]]
381			// [[/content]]
382			#endif
383