Uv sutherland

MSVStudio/thermopack.def

@@ -1,5 +1,5 @@
 ; File generated by thermopack/addon/pyUtils/exports/export_list.py
-; Time stamp: 2023-02-21T15:48:41.476484
+; Time stamp: 2023-02-27T13:17:37.477241
 ; @NAME@.def : Declares the module parameters.
 
 LIBRARY
@@ -73,6 +73,7 @@ EXPORTS
 	eoslibinit_mp_redefine_critical_parameters_
 	eoslibinit_mp_init_lj_
 	eoslibinit_mp_init_ljs_
+	eoslibinit_mp_init_uv_
 	eostv_mp_internal_energy_tv_
 	eostv_mp_entropy_tv_
 	eostv_mp_pressure_

Makefile.code

@@ -123,6 +123,8 @@ OBJE = \
        $(ODIR)/idealh2.o \
        $(ODIR)/ideal.o \
        $(ODIR)/cpa_parameters.o \
+       $(ODIR)/pair_potentials.o \
+       $(ODIR)/sutherland_a1a2a3.o \
        $(ODIR)/cbselect.o \
        $(ODIR)/hardsphere_bmcsl.o \
        $(ODIR)/pc_saft_nonassoc.o \
@@ -138,6 +140,7 @@ OBJE = \
        $(ODIR)/saftvrmie_interface.o \
        $(ODIR)/hardsphere_wca.o \
        $(ODIR)/lj_splined.o \
+       $(ODIR)/uv_theory.o \
        $(ODIR)/saft_globals.o \
        $(ODIR)/saft_rdf.o \
        $(ODIR)/saft_association.o \
@@ -384,6 +387,11 @@ $(ODIR)/cpa_parameters.o: $(SRC)/cpa_parameters.f90 \
 	$(SRC)/cubic_eos.f90 $(SRC)/assocschemeutils.f90 \
 	$(SRC)/stringmod.f90 $(SRC)/compdata.f90
 	@$(call make_object, $(SRC)/cpa_parameters.f90)
+$(ODIR)/pair_potentials.o: $(SRC)/pair_potentials.f90 \
+	$(SRC)/eosdata.f90 $(SRC)/thermopack_constants.f90 \
+	$(SRC)/hyperdual_mod.f90 $(SRC)/numconstants.f90 \
+	$(SRC)/nonlinear_solvers.f90
+	@$(call make_object, $(SRC)/pair_potentials.f90)
 $(ODIR)/critical.o: $(SRC)/critical.f90 \
 	$(SRC)/thermopack_constants.f90 $(SRC)/thermopack_var.f90 \
 	$(SRC)/eos.f90 $(SRC)/eostv.f90 \
@@ -415,6 +423,7 @@ $(ODIR)/eos_container.o: $(SRC)/eos_container.f90 \
 	$(SRC)/compdata.f90 $(SRC)/eos_parameters.f90 \
 	$(SRC)/cubic_eos.f90 $(SRC)/thermopack_var.f90 \
 	$(SRC)/saftvrmie_containers.f90 $(SRC)/lj_splined.f90 \
+	$(SRC)/uv_theory.f90 \
 	$(SRC)/pc_saft_nonassoc.f90 $(SRC)/extcsp.f90 \
 	$(SRC)/eosdata.f90 $(SRC)/pets.f90
 	@$(call make_object, $(SRC)/eos_container.f90)
@@ -478,6 +487,7 @@ $(ODIR)/external.o: $(SRC)/external.f90 \
 $(ODIR)/h2o_gibbs.o: $(SRC)/h2o_gibbs.f90
 	@$(call make_object, $(SRC)/h2o_gibbs.f90)
 $(ODIR)/hardsphere_bmcsl.o: $(SRC)/hardsphere_bmcsl.f90 \
+	$(SRC)/quadratures.f90 $(SRC)/hyperdual_mod.f90 \
 	$(SRC)/thermopack_constants.f90 $(SRC)/numconstants.f90
 	@$(call make_object, $(SRC)/hardsphere_bmcsl.f90)
 $(ODIR)/hardsphere_wca.o: $(SRC)/hardsphere_wca.f90 \
@@ -519,14 +529,27 @@ $(ODIR)/leekesler.o: $(SRC)/leekesler.f90 \
 	$(SRC)/numconstants.f90
 	@$(call make_object, $(SRC)/leekesler.f90)
 $(ODIR)/lj_splined.o: $(SRC)/lj_splined.f90 \
-	$(SRC)/saftvrmie_containers.f90 $(SRC)/saftvrmie_options.f90 \
-	$(SRC)/saftvrmie_utils.f90 $(SRC)/numconstants.f90 \
-	$(SRC)/eosdata.f90 $(SRC)/thermopack_constants.f90 \
-	$(SRC)/thermopack_var.f90 $(SRC)/hardsphere_wca.f90 \
-	$(SRC)/saftvrmie_dispersion.f90 $(SRC)/saftvrmie_hardsphere.f90 \
-	$(SRC)/compdata.f90 $(SRC)/stringmod.f90 \
-	$(SRC)/lj_splined.f90
+	$(SRC)/saftvrmie_containers.f90 $(SRC)/saftvrmie_utils.f90 \
+	$(SRC)/numconstants.f90 $(SRC)/eosdata.f90 \
+	$(SRC)/thermopack_constants.f90 $(SRC)/thermopack_var.f90 \
+	$(SRC)/hardsphere_wca.f90 $(SRC)/saftvrmie_dispersion.f90 \
+	$(SRC)/saftvrmie_hardsphere.f90 $(SRC)/compdata.f90 \
+	$(SRC)/stringmod.f90 $(SRC)/lj_splined.f90 \
+	$(SRC)/saft_interface.f90 $(SRC)/saftvrmie_options.f90
 	@$(call make_object, $(SRC)/lj_splined.f90)
+$(ODIR)/sutherland_a1a2a3.o: $(SRC)/sutherland_a1a2a3.f90 \
+	$(SRC)/hyperdual_mod.f90 $(SRC)/numconstants.f90 \
+	$(SRC)/thermopack_constants.f90	 
+	@$(call make_object, $(SRC)/sutherland_a1a2a3.f90)
+$(ODIR)/uv_theory.o: $(SRC)/uv_theory.f90 \
+	$(SRC)/saftvrmie_containers.f90 $(SRC)/saftvrmie_utils.f90 \
+	$(SRC)/numconstants.f90 $(SRC)/eosdata.f90 \
+	$(SRC)/thermopack_constants.f90 $(SRC)/thermopack_var.f90 \
+	$(SRC)/hardsphere_wca.f90 $(SRC)/saftvrmie_dispersion.f90 \
+	$(SRC)/compdata.f90 $(SRC)/pair_potentials.f90 \
+	$(SRC)/stringmod.f90 $(SRC)/saft_interface.f90 \
+	$(SRC)/saftvrmie_options.f90
+	@$(call make_object, $(SRC)/uv_theory.f90)
 $(ODIR)/mbwr.o: $(SRC)/mbwr.f90 \
 	$(SRC)/thermopack_constants.f90 $(SRC)/mbwrdata.f90 \
 	$(SRC)/compdatadb.f90 $(SRC)/compdata_init.f90 \

addon/pyExamples/uv-theory.py

@@ -0,0 +1,46 @@
+#!/usr/bin/python
+#Modify system path
+import sys
+sys.path.insert(0,'../pycThermopack/')
+# Importing pyThermopack
+from thermopack.uv_theory import uv_theory
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Instantiate and init uv_theory object. Model by van Westen and Gross (10.1063/5.0073572)
+
+def plot_BH_WCA_comparison(ax, comp, pert, paramref, kw):
+    uv = uv_theory(comp, pert, paramref)
+    uv.set_tmin(5.0)
+
+    if paramref == "SUTHERLAND":
+        print ("set params")
+        uv.set_params(ic=1,jc=1,nt=2, \
+                      c_vec=[3.8847749838093475, -3.8847749838093475], \
+                      lam_vec=[12.26, 6.0], \
+                      sigma=3.41e-10, \
+                      epsdivk=118.7, \
+                      beta_expo=None)
+
+    # Plot phase envelope
+    z = np.array([1.0])
+    T, P, v = uv.get_envelope_twophase(1.0e4, z, maximum_pressure=1.5e7, calc_v=True)
+    Tc, vc, Pc = uv.critical(z)
+    ax.plot(1/v, P*1.0e-6, **kw)
+    ax.plot([1/vc], [Pc*1.0e-6], "ko")
+
+ax = plt.gca()
+plot_BH_WCA_comparison(ax, "Ar", "WCA", "SVRMIE", kw={'label':"UV-WCA", 'ls':'-'})
+plot_BH_WCA_comparison(ax, "Ar", "BH", "SVRMIE", kw={'label':"UV-BH", 'ls':'--'})
+
+plot_BH_WCA_comparison(ax, "Ar", "BH", "SUTHERLAND", kw={'label':"UV-LAFITTE", 'ls':'-.'})
+
+ax.set_ylabel(r"$P$ (MPa)")
+ax.set_xlabel(r"$\rho$ (mol/m3)")
+ax.set_title("Argon UV-theory phase diagram")
+ax.set_xlim(70,35000)
+ax.set_ylim(0,6)
+ax.grid()
+plt.legend()
+plt.savefig("AR_COMPARISON.pdf")
+plt.show()

addon/pyUtils/exports/export_list.py

@@ -117,6 +117,13 @@ def append_export(method, module=None, isBindC=False):
 append_export("thermopack_setwsparam")
 append_export("thermopack_get_volume_shift_parameters")
 append_export("thermopack_set_volume_shift_parameters")
+append_export("thermopack_setclassicfitparam")
+append_export("thermopack_settwuparam")
+append_export("thermopack_gettwuparam")
+append_export("thermopack_set_alpha_corr")
+
+
+
 
 append_export("get_bp_term", "binaryplot")
 append_export("vllebinaryxy", "binaryplot")
@@ -156,6 +163,7 @@ def append_export(method, module=None, isBindC=False):
 append_export("redefine_critical_parameters", "eoslibinit")
 append_export("init_lj", "eoslibinit")
 append_export("init_ljs", "eoslibinit")
+append_export("init_uv", "eoslibinit")
 
 append_export("internal_energy_tv", "eostv")
 append_export("entropy_tv", "eostv")
@@ -201,6 +209,9 @@ def append_export(method, module=None, isBindC=False):
 
 append_export("lng_ii_pc_saft_tvn", "pc_saft_nonassoc")
 
+append_export("reset_sutsum_external_ij", "uv_theory")
+
+
 append_export("setphtolerance", "ph_solver")
 append_export("twophasephflash", "ph_solver")
 

addon/pycThermopack/thermopack/cpa.py

@@ -107,7 +107,7 @@ def get_kij(self, c1, c2):
             c2 (int): Component two
 
         Returns:
-            kij (array_like): i-j interaction parameter (2 parameters)
+            kij (array_like): i-j interaction parameter (a_kij, epskl_kij)
         """
         self.activate()
         c1_c = c_int(c1)
@@ -131,7 +131,7 @@ def set_kij(self, c1, c2, kij):
         Args:
             c1 (int): Component one
             c2 (int): Component two
-            kij (array_like): i-j interaction parameter (2 parameters)
+            kij (array_like): i-j interaction parameters (a_kij, epskl_kij)
         """
         self.activate()
         c1_c = c_int(c1)

addon/pycThermopack/thermopack/cubic.py

@@ -59,6 +59,8 @@ def __init__(self, comps=None, eos=None, mixing="vdW", alpha="Classic",
         self.s_get_covolumes = getattr(self.tp, self.get_export_name("cubic_eos", "get_covolumes"))
         self.s_get_energy_constants = getattr(self.tp, self.get_export_name("cubic_eos", "get_energy_constants"))
 
+        self.s_set_alpha_corr = getattr(self.tp, self.get_export_name("", "thermopack_set_alpha_corr"))
+
         if None not in (comps, eos):
             self.init(comps, eos, mixing, alpha, parameter_reference, volume_shift)
         elif self is cubic:
@@ -502,3 +504,32 @@ def get_energy_constants(self):
         self.s_get_energy_constants.restype = None
         self.s_get_energy_constants(a_c)
         return np.array(a_c)
+
+
+    def set_alpha_corr(self, ic, corrname, coeffs):
+        """Set alpha correlation
+
+        Args:
+            ic (in): Component number
+            corrname (string): Name of correlation
+            coeffs (ndarray): Coefficients in correlation
+        """
+
+        numparam_c = c_int(len(coeffs))
+        ic_c = c_int(ic)
+        corrname_string_c = c_char_p(corrname.strip().encode('ascii'))
+        corrname_string_len_c = c_len_type(len(corrname))
+        coeffs_c = (c_double * len(coeffs)) (*coeffs)
+        self.s_set_alpha_corr.argtypes = [POINTER(c_int),
+                                          POINTER(c_int),
+                                          c_char_p,
+                                          POINTER(c_double),
+                                          c_len_type]
+
+        self.s_set_alpha_corr.restype = None
+
+        self.s_set_alpha_corr(numparam_c,
+                              ic_c,
+                              corrname_string_c,
+                              coeffs_c,
+                              corrname_string_len_c)

addon/pycThermopack/thermopack/uv_theory.py

@@ -0,0 +1,117 @@
+from ctypes import *
+import numpy as np
+from . import thermo
+
+c_len_type = thermo.c_len_type
+
+
+class uv_theory(thermo.thermo):
+    """
+    Interface to UV-theory
+    """
+    def __init__(self, comps=None, model="BH", parameter_reference="Default"):
+        """
+        Initialize UV-Mie theory EoS
+        Args:
+            comps (str): Comma separated list of component names
+            model (str, optional): "BH" (default), "WCA"
+            parameter_reference (str, optional): Which parameters to use?. Defaults to "Default".
+        """
+        # Load dll/so
+        thermo.thermo.__init__(self)
+
+        # Init methods
+        self.s_eoslibinit_init_uv = getattr(self.tp, self.get_export_name("eoslibinit", "init_uv"))
+        self.s_set_sutsum_params = getattr(self.tp, self.get_export_name(
+            "uv_theory", "reset_sutsum_external_ij"))
+
+        if comps is not None:
+            self.init(comps, model, parameter_reference)
+
+    #################################
+    # Init
+    #################################
+    def init(self, comps, model="BH", parameter_reference="Default"):
+        """Initialize UV theory for Mie fluids
+
+        Args:
+            comps (str): Comma separated list of component names
+            model (str, optional): "BH" (default), "WCA"
+            parameter_reference (str, optional): Which parameters to use?. Defaults to "Default".
+        """
+        self.activate()
+        comp_string_c = c_char_p(comps.encode('ascii'))
+        comp_string_len = c_len_type(len(comps))
+        model_string_c = c_char_p(model.encode('ascii'))
+        model_string_len = c_len_type(len(model))
+        ref_string_c = c_char_p(parameter_reference.encode('ascii'))
+        ref_string_len = c_len_type(len(parameter_reference))
+
+        self.s_eoslibinit_init_uv.argtypes = [c_char_p,
+                                              c_char_p,
+                                              c_char_p,
+                                              c_len_type,
+                                              c_len_type,
+                                              c_len_type]
+
+        self.s_eoslibinit_init_uv.restype = None
+
+        self.s_eoslibinit_init_uv(comp_string_c,
+                                  model_string_c,
+                                  ref_string_c,
+                                  comp_string_len,
+                                  model_string_len,
+                                  ref_string_len)
+
+        self.nc = max(len(comps.split(" ")),len(comps.split(",")))
+
+
+    #def set_params(self, ic, m, sigma, eps_div_k, lambda_a, lambda_r):
+    def set_params(self,ic,jc,nt,c_vec,lam_vec,sigma,epsdivk,beta_expo=None):
+        # """Set pure fluid parameters
+
+        # Args:
+        #     ic (int): Component index
+        #     m (float): Mean number of segments.
+        #     sigma (float): Temperature-independent segment diameter [m].
+        #     eps_div_k (float): Well depth divided by Boltzmann's const. [K].
+        #     lambda_a (float): Attractive exponent of the Mie potential
+        #     lambda_r (float): Repulsive exponent of the Mie potential
+        # """
+        self.activate()
+
+        ic_c = c_int(ic)
+        jc_c = c_int(jc)
+        nt_c = c_int(nt)
+
+        cvec_c   = (c_double * nt)(*c_vec)
+        lamvec_c = (c_double * nt)(*lam_vec)
+
+        sigma_c = c_double(sigma)
+        eps_c = c_double(epsdivk)
+
+        null_pointer = POINTER(c_int)()
+        if beta_expo is None:
+            bexp_c = null_pointer
+        else:
+            bexp_c = (c_int * nt)(*beta_expo)
+
+        self.s_set_sutsum_params.argtypes = [POINTER(c_int),
+                                             POINTER(c_int),
+                                             POINTER(c_int),
+                                             POINTER(c_double),
+                                             POINTER(c_double),
+                                             POINTER(c_double),
+                                             POINTER(c_double),
+                                             POINTER(c_int)]
+
+        self.s_set_sutsum_params.restype = None
+
+        self.s_set_sutsum_params(byref(ic_c),
+                                 byref(jc_c),
+                                 byref(nt_c),
+                                 cvec_c,
+                                 lamvec_c,
+                                 byref(sigma_c),
+                                 byref(eps_c),
+                                 bexp_c)

libthermopack_export.symbols

@@ -1,5 +1,5 @@
 # File generated by thermopack/addon/pyUtils/exports/export_list.py
-# Time stamp: 2023-02-21T15:48:41.476345
+# Time stamp: 2023-03-28T10:48:56.952855
 # Explicit list of symbols to be exported
 	_thermopack_init_c
 	_thermopack_pressure_c
@@ -36,6 +36,10 @@
 	_thermopack_setwsparam_
 	_thermopack_get_volume_shift_parameters_
 	_thermopack_set_volume_shift_parameters_
+	_thermopack_setclassicfitparam_
+	_thermopack_settwuparam_
+	_thermopack_gettwuparam_
+	_thermopack_set_alpha_corr_
 	*get_bp_term*
 	*vllebinaryxy*
 	*global_binary_plot*
@@ -69,6 +73,7 @@
 	*redefine_critical_parameters*
 	*init_lj*
 	*init_ljs*
+	*init_uv*
 	*internal_energy_tv*
 	*entropy_tv*
 	*pressure*
@@ -103,6 +108,7 @@
 	*ljs_wca_model_control*
 	*ljs_wca_set_pure_params*
 	*ljs_wca_get_pure_params*
+	*reset_sutsum_external_ij*
 	*fres_multipol*
 	*multipol_model_control*
 	*lng_ii_pc_saft_tvn*