Implemented the Batzle and Wang 1992 viscosity EOS in the source code…

… + added new option in Diamond to activate it.

ICFERST/schemas/diagnostic_algorithms.rnc

@@ -1089,6 +1089,17 @@ tensor_python_diagnostic_algorithm =
          }?
       }
    )
+viscosity_EOS_algorithm =
+   (
+        ## Viscosity equation of state.
+        ## This diagnostic is internal - i.e. it is calculated somewhere within
+        ## the main code, and is not wrapped by the automatic diagnostics
+        ## wrappers.
+        element viscosity_BW {
+           attribute name { "Internal" },
+           attribute material_phase_support { "multiple" }
+        }
+     )
 
 # Adaptivity diagnostics
 scalar_edge_lengths_algorithm =
@@ -1707,4 +1718,5 @@ tensor_diagnostic_algorithms |= helmholtz_smoothed_tensor_algorithm
 tensor_diagnostic_algorithms |= helmholtz_anisotropic_smoothed_tensor_algorithm
 tensor_diagnostic_algorithms |= field_tolerance_algorithm
 tensor_diagnostic_algorithms |= tensor_python_diagnostic_algorithm
+tensor_diagnostic_algorithms |= viscosity_EOS_algorithm
 tensor_diagnostic_algorithms |= tensor_difference_algorithm

ICFERST/schemas/diagnostic_algorithms.rng

@@ -145,6 +145,27 @@ wrappers.</a:documentation>
       </attribute>
     </element>
   </define>
+  <define name="viscosity_EOS_algorithm">
+    <element name="viscosity_EOS">
+      <a:documentation>Viscosity equation of state.
+        This diagnostic is internal - i.e. it is calculated somewhere within
+  the main code, and is not wrapped by the automatic diagnostics
+  wrappers.</a:documentation>
+      <choice>
+        <element name="viscosity_BW">
+          <a:documentation>Batzle and Wang (1992) viscosity formulation.
+        Use for basin scale simulations.
+        Valid for ranges: T in [0,250] degrees Celsius, C in [0,0.46] kg/kg salinity.</a:documentation>
+          <attribute name="name">
+            <value>Internal</value>
+          </attribute>
+          <attribute name="material_phase_support">
+            <value>multiple</value>
+          </attribute>
+        </element>
+      </choice>
+    </element>
+  </define>
   <define name="legacy_internal_algorithm">
     <element name="algorithm">
       <a:documentation>This diagnostic is deprecated - i.e. it has been replaced by an
@@ -2409,6 +2430,9 @@ integral being zero.</a:documentation>
   <define name="tensor_diagnostic_algorithms" combine="choice">
     <ref name="tensor_python_diagnostic_algorithm"/>
   </define>
+  <define name="tensor_diagnostic_algorithms" combine="choice">
+    <ref name="viscosity_EOS_algorithm"/>
+  </define>
   <define name="tensor_diagnostic_algorithms" combine="choice">
     <ref name="tensor_difference_algorithm"/>
   </define>

ICFERST/schemas/multiphase.rnc

@@ -4092,6 +4092,7 @@ surface_integral_stats_vector.surface_integral &=
    )
 
 phase_viscosity =(
+                  ## Use this option to specify the viscosity of fluids
                   element tensor_field {
                      attribute name { "Viscosity" },
                      attribute rank { "2" },
@@ -4100,13 +4101,13 @@ phase_viscosity =(
                             pressure_mesh_choice,
                             prescribed_tensor_field_no_adapt
                         }|
-                        ## For electrical modelling only - holds electrical conductivity
                         element diagnostic {
                            (
                               tensor_python_diagnostic_algorithm |
-                              internal_algorithm
+                              internal_algorithm|
+                              viscosity_EOS_algorithm
                            ),
-                           velocity_mesh_choice,
+                           pressure_mesh_choice,
                            diagnostic_tensor_field
                         }
                      )
@@ -4252,7 +4253,7 @@ element Density{(
                         }|
                         ## Batzle and Wang (1992) EOS
                         ## Use for basin scale simulations.
-                        ## Valid for ranges: 5<=P<=100 MPa, 20<=T<=350 degrees Celsius, C<=0.32 kg/kg salinity.
+                        ## Valid for ranges: P in [5,100] MPa, T in [20,350] degrees Celsius, C in [0,0.32] kg/kg salinity.
                         ##
                         ## NOTE: Reference density is calculated using reference C0, T0, P0.
                         element BW_eos {
@@ -4302,7 +4303,7 @@ phase_properties_options =
 
                 ## Viscosity of the current phase (tensor field)
                 element Viscosity {
-                phase_viscosity,
+                    phase_viscosity,
                     ## Linearises a P2 field. Use this if having stability issues with high-order discretisations
                     element linearise_viscosity {empty}?,
                element viscosity_scheme {

ICFERST/schemas/multiphase.rng

@@ -4979,6 +4979,7 @@ Defaults to zero at machine tolerance (epsilon(0.0)) if not selected.</a:documen
   </define>
   <define name="phase_viscosity">
     <element name="tensor_field">
+      <a:documentation>Viscosity of the current phase (tensor field)</a:documentation>
       <attribute name="name">
         <value>Viscosity</value>
       </attribute>
@@ -4991,12 +4992,12 @@ Defaults to zero at machine tolerance (epsilon(0.0)) if not selected.</a:documen
           <ref name="prescribed_tensor_field_no_adapt"/>
         </element>
         <element name="diagnostic">
-          <a:documentation>For electrical modelling only - holds electrical conductivity</a:documentation>
           <choice>
             <ref name="tensor_python_diagnostic_algorithm"/>
             <ref name="internal_algorithm"/>
+            <ref name="viscosity_EOS_algorithm"/>
           </choice>
-          <ref name="velocity_mesh_choice"/>
+          <ref name="pressure_mesh_choice"/>
           <ref name="diagnostic_tensor_field"/>
         </element>
       </choice>
@@ -5167,11 +5168,11 @@ Introduce the same name as the scalar field being used. For example: Tracer_CO2<
               </zeroOrMore>
             </element>
             <element name="BW_eos">
-              <a:documentation>Batzle and Wang (1992) EOS
-Use for basin scale simulations.
-Valid for ranges: 5&lt;=P&lt;=100 MPa, 20&lt;=T&lt;=350 degrees Celsius, C&lt;=0.32 kg/kg salinity.
+              <a:documentation>Batzle and Wang (1992) EOS.
+  Use for basin scale simulations.
+  Valid for ranges: P in [5,100] MPa, T in [20,350] degrees Celsius, C in [0,0.32] kg/kg salinity.
 
-NOTE: Reference density is calculated using reference C0, T0, P0.</a:documentation>
+  NOTE: Reference density is calculated using reference C0, T0, P0.</a:documentation>
               <optional>
                 <element name="C0">
                   <a:documentation>Reference concentration (default 0)</a:documentation>
@@ -5235,8 +5236,8 @@ Component properties overwrite the ones given for a phase and are mandatory for
         <a:documentation>Phase density</a:documentation>
       </ref>
       <element name="Viscosity">
-        <a:documentation>Viscosity of the current phase (tensor field)</a:documentation>
-        <ref name="phase_viscosity"/>
+        <a:documentation>Use this option to specify the viscosity of fluids</a:documentation>
+          <ref name="phase_viscosity"/>
         <optional>
           <element name="linearise_viscosity">
             <a:documentation>Linearises a P2 field. Use this if having stability issues with high-order discretisations</a:documentation>

ICFERST/src/Multiphase_TimeLoop.F90

@@ -250,6 +250,7 @@ subroutine MultiFluids_SolveTimeLoop( state, &
         integer :: phreeqc_id
         double precision, ALLOCATABLE, dimension(:,:) :: concetration_phreeqc
         real :: total_mass_metal_before_adapt, total_mass_metal_after_adapt, total_mass_metal_after_correction, total_mass_metal_after_bound
+        logical :: viscosity_EOS
 #ifdef HAVE_ZOLTAN
       real(zoltan_float) :: ver
       integer(zoltan_int) :: ierr
@@ -439,6 +440,8 @@ subroutine MultiFluids_SolveTimeLoop( state, &
             end if
         end do
 
+        call compute_viscosity_EOS( state, Mdims )
+
         if( have_option( '/material_phase[0]/multiphase_properties/capillary_pressure' ) ) &
             have_extra_DiffusionLikeTerm = .true.
         if ( have_option( '/mesh_adaptivity/hr_adaptivity' ) ) then
@@ -909,6 +912,8 @@ subroutine MultiFluids_SolveTimeLoop( state, &
             !Time to compute the self-potential if required
             if (write_all_stats .and. have_option("/porous_media/SelfPotential")) &
                     call Assemble_and_solve_SP(Mdims, state, packed_state, ndgln, Mmat, Mspars, CV_funs, CV_GIdims)
+            ! Compute viscosity
+            call compute_viscosity_EOS( state, Mdims )
             !Now compute diagnostics
             call calculate_diagnostic_variables( state, exclude_nonrecalculated = .true. )
             !calculate_diagnostic_variables_new <= computes other diagnostics such as python-based fields

ICFERST/src/multi_eos.F90

@@ -872,7 +872,7 @@ subroutine Calculate_PorousMedia_AbsorptionTerms( nphase, state, packed_state, P
        real :: Angle, bad_element_perm_mult, Max_aspect_ratio, height ! the height of an isosceles triangle for the top angle to be equal to the trigger angle
        real, dimension(Mdims%ndim,Mdims%ndim) :: trans_matrix, rot_trans_matrix ! for bad_element permeability transformation matrix
        real, parameter :: pi = acos(0.0) * 2.0 ! Define pi
-       character( len = option_path_len ) :: option_path_python
+       character( len = option_path_len ) :: option_path_python, option_path_viscosity_EOS
 
        perm => extract_tensor_field( packed_state, "Permeability" )
        !Define n_in_pres based on the local version of nphase
@@ -888,13 +888,21 @@ subroutine Calculate_PorousMedia_AbsorptionTerms( nphase, state, packed_state, P
        DO IPHASE = 1, Mdims%nphase!Get viscosity for all the phases
         option_path_python = "/material_phase["// int2str( iphase - 1 )//"]/phase_properties/Viscosity/tensor_field"//&
         "::Viscosity/diagnostic/algorithm::tensor_python_diagnostic"
+        option_path_viscosity_EOS = "/material_phase["// int2str( iphase - 1 )//"]/phase_properties/Viscosity/tensor_field"//&
+        "::Viscosity/diagnostic/viscosity_EOS"
           if (have_option(trim(option_path_python)))then
             state_viscosity => extract_tensor_field( state( iphase ), 'Viscosity' )
             call multi_compute_python_field(state, iphase, trim( option_path_python ), tfield = state_viscosity)
             !Copy into state
             do i = 1, Mdims%cv_nonods
               viscosities(iphase,i) = state_viscosity%val(1,1,i)
             end do
+          else if (have_option(trim(option_path_viscosity_EOS))) then
+            call compute_viscosity_EOS( state, Mdims )
+            state_viscosity => extract_tensor_field( state( iphase ), 'Viscosity' )
+            do i = 1, Mdims%cv_nonods
+              viscosities(iphase,i) = state_viscosity%val(1,1,i)
+            end do
           else
             call set_viscosity(nphase, Mdims, state, viscosities)
           end if
@@ -1944,6 +1952,75 @@ subroutine calculate_viscosity( state, Mdims, ndgln, Momentum_Diffusion, Momentu
 
     end subroutine calculate_viscosity
 
+    subroutine compute_viscosity_EOS( state, Mdims )
+      implicit none
+      type( multi_dimensions ), intent( in ) :: Mdims
+      type( state_type ), dimension( : ), intent( inout ) :: state
+      type( tensor_field ), pointer :: t_field
+      integer :: iphase, stat, cv_nod
+      type( scalar_field ), pointer :: temperature, concentration
+      logical :: viscosity_BW, have_temperature_field, have_concentration_field
+
+        do iphase = 1, Mdims%nphase
+          viscosity_BW = have_option("/material_phase["// int2str( iphase - 1 )//"]/phase_properties/Viscosity/tensor_field"//&
+          "::Viscosity/diagnostic/viscosity_EOS/viscosity_BW::Internal")
+          if (viscosity_BW) then
+            temperature => extract_scalar_field( state( iphase ), 'Temperature', stat )
+            have_temperature_field = ( stat == 0 )
+            Concentration => extract_scalar_field( state( iphase ), 'Concentration', stat )
+            have_concentration_field = ( stat == 0 )
+            t_field => extract_tensor_field( state( iphase ), 'Viscosity', stat )
+            if (.not. (have_temperature_field)) then
+              FLAbort( "Temperature field needed for BW1992 viscosity EOS." )
+            else
+              do cv_nod=1,Mdims%cv_nonods
+                if (temperature%val(cv_nod) < 273.15) then
+                  t_field%val(1, 1, cv_nod) = 1.e-3
+                  t_field%val(1, 2, cv_nod) = 0.0
+                  t_field%val(1, 3, cv_nod) = 0.0
+                  t_field%val(2, 1, cv_nod) = 0.0
+                  t_field%val(2, 2, cv_nod) = 1.e-3
+                  t_field%val(2, 3, cv_nod) = 0.0
+                  t_field%val(3, 1, cv_nod) = 0.0
+                  t_field%val(3, 2, cv_nod) = 0.0
+                  t_field%val(3, 3, cv_nod) = 1.e-3
+                else
+                  if (have_concentration_field) then
+                    t_field%val(1, 1, cv_nod) = 1e-3 * (0.1 + 0.333*concentration%val(cv_nod) + (1.65 + 91.9 * (concentration%val(cv_nod))**3) * exp(-(0.42 * ((concentration%val(cv_nod))**0.8-0.17)**2 + 0.045) &
+                    * (temperature%val(cv_nod) - 273.15)**0.8))
+                    t_field%val(1, 2, cv_nod) = 0.0
+                    t_field%val(1, 3, cv_nod) = 0.0
+                    t_field%val(2, 1, cv_nod) = 0.0
+                    t_field%val(2, 2, cv_nod) = 1e-3 * (0.1 + 0.333*concentration%val(cv_nod) + (1.65 + 91.9 * (concentration%val(cv_nod))**3) * exp(-(0.42 * ((concentration%val(cv_nod))**0.8-0.17)**2 + 0.045) &
+                    * (temperature%val(cv_nod) - 273.15)**0.8))
+                    t_field%val(2, 3, cv_nod) = 0.0
+                    t_field%val(3, 1, cv_nod) = 0.0
+                    t_field%val(3, 2, cv_nod) = 0.0
+                    t_field%val(3, 3, cv_nod) = 1e-3 * (0.1 + 0.333*concentration%val(cv_nod) + (1.65 + 91.9 * (concentration%val(cv_nod))**3) * exp(-(0.42 * ((concentration%val(cv_nod))**0.8-0.17)**2 + 0.045) &
+                    * (temperature%val(cv_nod) - 273.15)**0.8))
+                  else
+                    t_field%val(1, 1, cv_nod) = 1e-3 * (0.1 + (1.65) * exp(-(0.42 * (-0.17)**2 + 0.045) * (temperature%val(cv_nod) - 273.15)**0.8))
+                    t_field%val(1, 2, cv_nod) = 0.0
+                    t_field%val(1, 3, cv_nod) = 0.0
+                    t_field%val(2, 1, cv_nod) = 0.0
+                    t_field%val(2, 2, cv_nod) = 1e-3 * (0.1 + (1.65) * exp(-(0.42 * (-0.17)**2 + 0.045) * (temperature%val(cv_nod) - 273.15)**0.8))
+                    t_field%val(2, 3, cv_nod) = 0.0
+                    t_field%val(3, 1, cv_nod) = 0.0
+                    t_field%val(3, 2, cv_nod) = 0.0
+                    t_field%val(3, 3, cv_nod) = 1e-3 * (0.1 + (1.65) * exp(-(0.42 * (-0.17)**2 + 0.045) * (temperature%val(cv_nod) - 273.15)**0.8))
+                  end if
+                end if
+              end do
+            end if
+
+            ! Make sure viscosity stays between bounds.
+            t_field%val = max(min(t_field%val,1.e-3), 1.e-4)
+          end if
+        end do
+
+    end subroutine compute_viscosity_EOS
+
+
 
 
     !sprint_to_do, re-use material_absoprtion by updating the values of the input absoprtion