From a786d0b62a09ca77f9a532aee8a93cdce55261be Mon Sep 17 00:00:00 2001 From: Meissam Bahlali Date: Wed, 17 Jul 2024 13:56:33 +0100 Subject: [PATCH] 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 | 12 ++++ ICFERST/schemas/diagnostic_algorithms.rng | 24 +++++++ ICFERST/schemas/multiphase.rnc | 11 ++-- ICFERST/schemas/multiphase.rng | 17 ++--- ICFERST/src/Multiphase_TimeLoop.F90 | 5 ++ ICFERST/src/multi_eos.F90 | 79 ++++++++++++++++++++++- 6 files changed, 134 insertions(+), 14 deletions(-) diff --git a/ICFERST/schemas/diagnostic_algorithms.rnc b/ICFERST/schemas/diagnostic_algorithms.rnc index 2002bd3c5..af3707705 100644 --- a/ICFERST/schemas/diagnostic_algorithms.rnc +++ b/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 diff --git a/ICFERST/schemas/diagnostic_algorithms.rng b/ICFERST/schemas/diagnostic_algorithms.rng index 573ec9224..768c9870f 100644 --- a/ICFERST/schemas/diagnostic_algorithms.rng +++ b/ICFERST/schemas/diagnostic_algorithms.rng @@ -145,6 +145,27 @@ wrappers. + + + 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. + + + 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. + + Internal + + + multiple + + + + + This diagnostic is deprecated - i.e. it has been replaced by an @@ -2409,6 +2430,9 @@ integral being zero. + + + diff --git a/ICFERST/schemas/multiphase.rnc b/ICFERST/schemas/multiphase.rnc index 1aedf6e00..c0878ddff 100755 --- a/ICFERST/schemas/multiphase.rnc +++ b/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 { diff --git a/ICFERST/schemas/multiphase.rng b/ICFERST/schemas/multiphase.rng index 3ade03939..362d385c7 100755 --- a/ICFERST/schemas/multiphase.rng +++ b/ICFERST/schemas/multiphase.rng @@ -4979,6 +4979,7 @@ Defaults to zero at machine tolerance (epsilon(0.0)) if not selected. + Viscosity of the current phase (tensor field) Viscosity @@ -4991,12 +4992,12 @@ Defaults to zero at machine tolerance (epsilon(0.0)) if not selected. - For electrical modelling only - holds electrical conductivity + - + @@ -5167,11 +5168,11 @@ Introduce the same name as the scalar field being used. For example: Tracer_CO2< - 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. + 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. + NOTE: Reference density is calculated using reference C0, T0, P0. Reference concentration (default 0) @@ -5235,8 +5236,8 @@ Component properties overwrite the ones given for a phase and are mandatory for Phase density - Viscosity of the current phase (tensor field) - + Use this option to specify the viscosity of fluids + Linearises a P2 field. Use this if having stability issues with high-order discretisations diff --git a/ICFERST/src/Multiphase_TimeLoop.F90 b/ICFERST/src/Multiphase_TimeLoop.F90 index 4168e59cd..c26bf6e72 100644 --- a/ICFERST/src/Multiphase_TimeLoop.F90 +++ b/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 diff --git a/ICFERST/src/multi_eos.F90 b/ICFERST/src/multi_eos.F90 index 5ccce0134..5f3307858 100644 --- a/ICFERST/src/multi_eos.F90 +++ b/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,6 +888,8 @@ 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) @@ -895,6 +897,12 @@ subroutine Calculate_PorousMedia_AbsorptionTerms( nphase, state, packed_state, P 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