Initial hydrate support

MSVStudio/thermopack.def

@@ -1,5 +1,5 @@
 ; File generated by thermopack/addon/pyUtils/exports/export_list.py
-; Time stamp: 2024-10-04T15:56:57.529901
+; Time stamp: 2024-12-01T11:39:41.420875
 ; @[email protected] : Declares the module parameters.
 
 LIBRARY
@@ -98,6 +98,10 @@ EXPORTS
 	eostv_mp_enthalpy_tvp_
 	fundamental_measure_theory_mp_fmt_energy_density_
 	hardsphere_bmcsl_mp_calc_bmcsl_gij_fmt_
+	init_hydrate_model_mp_hydrate_
+	fugacity_water_in_hydrate_tpx_mp_hydrate_
+	fugacity_water_in_hydrate_tvn_mp_hydrate_
+	map_hydrate_appearance_curve_mp_hydrate_curves_
 	ideal_mp_set_standard_entropy_
 	ideal_mp_get_standard_entropy_
 	ideal_mp_set_enthalpy_of_formation_
@@ -117,6 +121,7 @@ EXPORTS
 	lj_splined_mp_ljs_wca_model_control_
 	lj_splined_mp_ljs_wca_set_pure_params_
 	lj_splined_mp_ljs_wca_get_pure_params_
+	multi_phase_envelope_tv_mp_multi_phase_envelope_plot_tv_
 	mut_solver_mp_solve_mu_t_
 	mut_solver_mp_solve_lnf_t_
 	mut_solver_mp_map_meta_isotherm_
@@ -196,6 +201,7 @@ EXPORTS
 	saturation_curve_mp_envelopeplot_
 	saturation_curve_mp_envelope_isentrope_cross_
 	saturation_curve_mp_pure_fluid_saturation_wrapper_
+	saturation_tv_mp_envelope_plot_tv_
 	saturation_point_locators_mp_locate_saturation_property_
 	saturation_point_locators_mp_property_index_from_string_
 	saturation_point_locators_mp_sat_points_based_on_prop_
@@ -229,4 +235,4 @@ EXPORTS
 	thermopack_var_mp_get_eos_identification_
 	tp_solver_mp_twophasetpflash_
 	sv_solver_mp_twophasesvflash_
-	uv_solver_mp_twophaseuvflash_
+	uv_solver_mp_twophaseuvflash_

Makefile.code

@@ -162,14 +162,16 @@ OBJE = \
        $(ODIR)/puresaturation.o \
        $(ODIR)/thermo_utils.o \
        $(ODIR)/critical.o \
+       $(ODIR)/stability.o \
        $(ODIR)/saturation.o \
        $(ODIR)/solideos.o \
        $(ODIR)/vls.o \
-       $(ODIR)/stability.o \
        $(ODIR)/tuning.o \
        $(ODIR)/tp_solver.o \
        $(ODIR)/state_functions.o \
        $(ODIR)/saturation_curve.o \
+       $(ODIR)/saturation_tv.o \
+       $(ODIR)/multi_phase_envelope_tv.o \
        $(ODIR)/solid_saturation.o \
        $(ODIR)/saturation_point_locators.o \
        $(ODIR)/joule_thompson_inversion.o \
@@ -183,6 +185,8 @@ OBJE = \
        $(ODIR)/binaryplot.o \
        $(ODIR)/speed_of_sound.o \
        $(ODIR)/saftvrmie_testing.o \
+       $(ODIR)/hydrate.o \
+       $(ODIR)/hydrate_curves.o \
        $(ODIR)/consistency.o \
        $(ODIR)/eoslibinit.o \
        $(ODIR)/complexmodelinit.o \
@@ -503,6 +507,23 @@ $(ODIR)/hardsphere_wca.o: $(SRC)/hardsphere_wca.f90 \
 	$(SRC)/thermopack_constants.f90 $(SRC)/hardsphere_wca.f90 \
 	$(SRC)/saftvrmie_options.f90
 	@$(call make_object, $(SRC)/hardsphere_wca.f90)
+$(ODIR)/hydrate.o: $(SRC)/hydrate.f90 \
+	$(SRC)/thermopack_constants.f90 $(SRC)/thermopack_var.f90 \
+	$(SRC)/eostv.f90 $(SRC)/saturation.f90 \
+	$(SRC)/stringmod.f90 $(SRC)/solideos.f90 \
+	$(SRC)/nonlinear_solvers.f90 $(SRC)/numconstants.f90 \
+	$(SRC)/quadratures.f90 $(SRC)/eos.f90
+	@$(call make_object, $(SRC)/hydrate.f90)
+$(ODIR)/hydrate_curves.o: $(SRC)/hydrate_curves.f90 \
+	$(SRC)/hydrate.f90 $(SRC)/thermopack_constants.f90 \
+	$(SRC)/thermopack_var.f90 $(SRC)/eostv.f90 \
+	$(SRC)/saturation.f90 $(SRC)/multi_phase_envelope_tv.f90 \
+	$(SRC)/saturation_tv.f90 $(SRC)/nonlinear_solvers.f90 \
+	$(SRC)/stability.f90 $(SRC)/numconstants.f90 \
+	$(SRC)/eosdata.f90 $(SRC)/cubic_eos.f90 \
+	$(SRC)/utilities.f90 $(SRC)/eos.f90 \
+	$(SRC)/thermo_utils.f90
+	@$(call make_object, $(SRC)/hydrate_curves.f90)
 $(ODIR)/hyperdual_mod.o: $(SRC)/hyperdual_mod.f90
 	@$(call make_object, $(SRC)/hyperdual_mod.f90)
 $(ODIR)/hyperdual_utility.o: $(SRC)/hyperdual_utility.f90 \
@@ -562,6 +583,16 @@ $(ODIR)/mbwrdata.o: $(SRC)/mbwrdata.f90 \
 $(ODIR)/mixdatadb.o: $(SRC)/mixdatadb.f90 \
 	$(SRC)/cubic_eos.f90 $(SRC)/assocschemeutils.f90
 	@$(call make_object, $(SRC)/mixdatadb.f90)
+$(ODIR)/multi_phase_envelope_tv.o: $(SRC)/multi_phase_envelope_tv.f90 \
+	$(SRC)/eos.f90 $(SRC)/eostv.f90 \
+	$(SRC)/thermopack_constants.f90 $(SRC)/nonlinear_solvers.f90 \
+	$(SRC)/thermopack_var.f90 $(SRC)/numconstants.f90 \
+	$(SRC)/thermo_utils.f90 $(SRC)/saturation.f90 \
+	$(SRC)/saturation_tv.f90 $(SRC)/saturation_curve.f90 \
+	$(SRC)/cubic_eos.f90 $(SRC)/utilities.f90 \
+	$(SRC)/eosdata.f90 $(SRC)/stability.f90 \
+	$(SRC)/critical.f90
+	@$(call make_object, $(SRC)/multi_phase_envelope_tv.f90)
 $(ODIR)/multiparameter_base.o: $(SRC)/multiparameter_base.f90 \
 	$(SRC)/numconstants.f90 $(SRC)/thermopack_constants.f90 \
 	$(SRC)/hyperdual_mod.f90
@@ -765,16 +796,17 @@ $(ODIR)/saturation.o: $(SRC)/saturation.f90 \
 	$(SRC)/eos.f90 $(SRC)/thermopack_constants.f90 \
 	$(SRC)/thermopack_var.f90 $(SRC)/nonlinear_solvers.f90 \
 	$(SRC)/numconstants.f90 $(SRC)/puresaturation.f90 \
-	$(SRC)/thermo_utils.f90 $(SRC)/eostv.f90
+	$(SRC)/thermo_utils.f90 $(SRC)/eostv.f90 \
+	$(SRC)/stability.f90
 	@$(call make_object, $(SRC)/saturation.f90)
 $(ODIR)/saturation_curve.o: $(SRC)/saturation_curve.f90 \
 	$(SRC)/eos.f90 $(SRC)/thermopack_constants.f90 \
 	$(SRC)/thermopack_var.f90 $(SRC)/numconstants.f90 \
 	$(SRC)/puresaturation.f90 $(SRC)/saturation.f90 \
 	$(SRC)/nonlinear_solvers.f90 $(SRC)/critical.f90 \
 	$(SRC)/thermo_utils.f90 $(SRC)/eostv.f90 \
-	$(SRC)/eosdata.f90 $(SRC)/utilities.f90 \
-	$(SRC)/solideos.f90
+	$(SRC)/eosdata.f90 $(SRC)/solideos.f90 \
+	$(SRC)/stability.f90 $(SRC)/utilities.f90
 	@$(call make_object, $(SRC)/saturation_curve.f90)
 $(ODIR)/saturation_point_locators.o: $(SRC)/saturation_point_locators.f90 \
 	$(SRC)/saturation.f90 $(SRC)/saturation_curve.f90 \
@@ -786,6 +818,16 @@ $(ODIR)/saturation_point_locators.o: $(SRC)/saturation_point_locators.f90 \
 	$(SRC)/thermo_utils.f90 $(SRC)/critical.f90 \
 	$(SRC)/solideos.f90 $(SRC)/vls.f90
 	@$(call make_object, $(SRC)/saturation_point_locators.f90)
+$(ODIR)/saturation_tv.o: $(SRC)/saturation_tv.f90 \
+	$(SRC)/eostv.f90 $(SRC)/thermopack_constants.f90 \
+	$(SRC)/thermopack_var.f90 $(SRC)/nonlinear_solvers.f90 \
+	$(SRC)/numconstants.f90 $(SRC)/saturation.f90 \
+	$(SRC)/saturation_curve.f90 $(SRC)/puresaturation.f90 \
+	$(SRC)/cubic_eos.f90 $(SRC)/utilities.f90 \
+	$(SRC)/thermo_utils.f90 $(SRC)/eos.f90 \
+	$(SRC)/eosdata.f90 $(SRC)/critical.f90 \
+	$(SRC)/stability.f90
+	@$(call make_object, $(SRC)/saturation_tv.f90)
 $(ODIR)/single_component.o: $(SRC)/single_component.f90 \
 	$(SRC)/eos_parameters.f90 $(SRC)/eosdata.f90 \
 	$(SRC)/mbwr_additional.f90 $(SRC)/compdata.f90 \

addon/pyExamples/hydrate_curves.py

@@ -0,0 +1,90 @@
+#!/usr/bin/python
+#Modify system path
+import sys
+sys.path.insert(0,'../pycThermopack/')
+from thermopack.cubic import cubic
+import numpy as np
+import matplotlib.pyplot as plt
+import itertools as it
+import csv
+
+cb = cubic("CO2,H2O","SRK")
+cb.init_hydrate()
+
+Z_H2O = np.array([30.0, 100.0, 250.0])*1e-6
+initial_pressure = 1.0e3
+minimum_temperature = 100.0
+maximum_pressure = 150.0e5
+p_scaling = 1.0e-6
+colors = [ "black", "blue", "red", "green", "m"]
+linestyles = [ "-", "--", ":", "-.", "."]
+
+z = np.zeros((2))
+for i in range(len(Z_H2O)):
+    z[0] = (1-Z_H2O[i])
+    z[1] = Z_H2O[i]
+    fluid, water = cb.get_multi_phase_envelope_tv(initial_pressure, z,
+                                                  minimum_temperature,
+                                                  maximum_pressure)
+    t_vals, p_vals = (fluid.t, fluid.p)
+    tw_vals, pw_vals = (water.t, water.p)
+    t_hyd_vals, p_hyd_vals = cb.get_hydrate_apperance_curve(initial_pressure, z,
+                                                            minimum_temperature,
+                                                            maximum_pressure)
+    plt.plot(t_vals, p_vals*p_scaling, linestyle="-", color=colors[i])
+    #plt.plot(tw_vals, pw_vals*p_scaling, linestyle=":", color=colors[i],
+    #         label="H2O: {} ppm".format(round(z[-1]*1e6)))
+    plt.plot(t_hyd_vals, p_hyd_vals*p_scaling, linestyle="--", color=colors[i],
+             label="Hyd. H2O: {} ppm".format(round(z[-1]*1e6)))
+
+# z = np.ones((2))*0.5
+# t_hyd_vals, p_hyd_vals = cb.get_hydrate_apperance_curve(initial_pressure, z,
+#                                                         minimum_temperature,
+#                                                         maximum_pressure)
+# plt.plot(t_hyd_vals, p_hyd_vals*p_scaling, linestyle="--", color="orange",
+#          label="Hyd. excess H2O")
+
+plt.title("SRK: CO2 and variable H2O")
+leg = plt.legend(loc="best", numpoints=1)
+leg.get_frame().set_linewidth(0.0)
+plt.ylabel(r"$P$ (MPa)")
+plt.xlabel(r"$T$ (K)")
+plt.tight_layout()
+plt.savefig("hydrate_co2.pdf")
+
+plt.figure()
+cb.init("CO2,N2,H2O","SRK","Classic","Classic")
+cb.init_hydrate()
+Z_CO2 = np.array([0.85,0.15])
+Z_H2O = np.array([30.0, 100.0, 250.0])*1e-6
+initial_pressure = 5.0e3
+
+z = np.zeros((3))
+for i in range(len(Z_H2O)):
+    z[0:2] = (1-Z_H2O[i])*Z_CO2
+    z[2] = Z_H2O[i]
+    fluid, water = cb.get_multi_phase_envelope_tv(initial_pressure, z,
+                                                  minimum_temperature,
+                                                  maximum_pressure)
+    t_vals, p_vals = (fluid.t, fluid.p)
+    tw_vals, pw_vals = (water.t, water.p)
+    t_hyd_vals, p_hyd_vals = cb.get_hydrate_apperance_curve(initial_pressure, z,
+                                                            minimum_temperature,
+                                                            maximum_pressure)
+
+    plt.plot(t_vals, p_vals*p_scaling, linestyle="-", color=colors[i])
+    # plt.plot(tw_vals, pw_vals*p_scaling, linestyle=":", color=colors[i],
+    #         label="H2O: {} ppm".format(round(z[-1]*1e6)))
+    plt.plot(t_hyd_vals, p_hyd_vals*p_scaling, linestyle="--", color=colors[i],
+             label="Hyd. H2O: {} ppm".format(round(z[-1]*1e6)))
+
+
+plt.title("SRK: CO2 (0.85), N2 (0.15) and variable H2O")
+leg = plt.legend(loc="best", numpoints=1)
+leg.get_frame().set_linewidth(0.0)
+plt.ylabel(r"$P$ (MPa)")
+plt.xlabel(r"$T$ (K)")
+plt.tight_layout()
+plt.savefig("hydrate_co2_n2.pdf")
+plt.show()
+plt.clf()

addon/pyExamples/multi_phase_curves.py

@@ -0,0 +1,88 @@
+#!/usr/bin/python
+#Modify system path
+import sys
+sys.path.insert(0,'../pycThermopack/')
+from thermopack.cubic import cubic
+import numpy as np
+import matplotlib.pyplot as plt
+import itertools as it
+import csv
+
+cb = cubic("CO2,N2,H2O","SRK")
+Z_CO2 = np.array([0.85,0.15])
+Z_H2O = np.array([0.05,1e-2,1e-3,1e-4])
+Z_H2O = np.array([5.0, 30.0, 100.0, 1000.0])*1e-6
+initial_pressure = 5.0e3
+minimum_temperature = 140.0
+maximum_pressure = 150.0e5
+p_scaling = 1.0e-6
+colors = [ "black", "blue", "red", "green"]
+linestyles = [ "-", "--", ":", "-."]
+
+z = np.zeros((3))
+for i in range(len(Z_H2O)):
+    z[0:2] = (1-Z_H2O[i])*Z_CO2
+    z[2] = Z_H2O[i]
+    fluid, water = cb.get_multi_phase_envelope_tv(initial_pressure, z,
+                                                  minimum_temperature,
+                                                  maximum_pressure)
+    t_vals, p_vals = (fluid.t, fluid.p)
+    tw_vals, pw_vals = (water.t, water.p)
+    plt.plot(t_vals, p_vals*p_scaling, linestyle="-", color=colors[i])
+    plt.plot(tw_vals, pw_vals*p_scaling, linestyle="--", color=colors[i],
+             label="H2O: {} ppm".format(round(z[-1]*1e6)))
+
+cb.init("CO2,N2","SRK","Classic","Classic")
+Z_CO2 = np.array([0.85,0.15])
+T, P, v = cb.get_envelope_twophase_tv(initial_pressure, Z_CO2,
+                                      maximum_pressure=maximum_pressure,
+                                      minimum_temperature=minimum_temperature)
+plt.plot(T, P*p_scaling, linestyle="-", color="cyan", label="No water")
+
+plt.title("SRK: CO2 (0.85), N2 (0.15) and variable H2O")
+leg = plt.legend(loc="best", numpoints=1)
+leg.get_frame().set_linewidth(0.0)
+plt.ylabel(r"$P$ (MPa)")
+plt.xlabel(r"$T$ (K)")
+
+plt.figure()
+cb.init("CO2,H2O","SRK","Classic","Classic")
+initial_pressure = 1.0e3
+minimum_temperature = 100.0
+maximum_pressure = 150.0e5
+p_scaling = 1.0e-6
+colors = [ "black", "blue", "red", "green"]
+linestyles = [ "-", "--", ":", "-."]
+
+arrays = []
+z = np.zeros((2))
+for i in range(len(Z_H2O)):
+    z[0] = (1-Z_H2O[i])
+    z[1] = Z_H2O[i]
+    fluid, water = cb.get_multi_phase_envelope_tv(initial_pressure, z,
+                                                  minimum_temperature,
+                                                  maximum_pressure)
+    t_vals, p_vals = (fluid.t, fluid.p)
+    tw_vals, pw_vals = (water.t, water.p)
+    arrays += [t_vals, p_vals, tw_vals, pw_vals]
+    plt.plot(t_vals, p_vals*p_scaling, linestyle="-", color=colors[i])
+    plt.plot(tw_vals, pw_vals*p_scaling, linestyle="--", color=colors[i],
+             label="H2O: {} ppm".format(round(z[-1]*1e6)))
+
+# with open('co2_h2o.csv', 'w') as f:
+#     csv.writer(f, delimiter='\t').writerows(it.zip_longest(*arrays, fillvalue =np.NaN))
+
+cb.init("CO2","SRK","Classic","Classic")
+Z_CO2 = np.array([1.0])
+T, P, v = cb.get_envelope_twophase_tv(initial_pressure, Z_CO2,
+                                       maximum_pressure=maximum_pressure,
+                                       minimum_temperature=minimum_temperature)
+plt.plot(T, P*p_scaling, linestyle="-", color="cyan", label="No water")
+
+plt.title("SRK: CO2 and variable H2O")
+leg = plt.legend(loc="best", numpoints=1)
+leg.get_frame().set_linewidth(0.0)
+plt.ylabel(r"$P$ (MPa)")
+plt.xlabel(r"$T$ (K)")
+
+plt.show()

addon/pyExamples/properties.py

@@ -46,7 +46,11 @@
 print("Compressibility factor: {}".format(zFac))
 
 sos = srk.speed_of_sound(T, P, x=z, y=z, z=z, betaV=0.0, betaL=1.0, phase=srk.LIQPH)
+sos_z = srk.speed_of_sound_tv(T, v, n=z)
+sos_n = srk.speed_of_sound_tv(T, v*2, n=z*2)
 print("Saturated liquid speed of sound: {} (m/s)".format(sos))
+print("Saturated liquid speed of sound: {} (m/s)".format(sos_z))
+print("Saturated liquid speed of sound: {} (m/s)".format(sos_n))
 
 # Properties evaluated in temperature volume and mol numbers
 p, dp  = srk.pressure_tv(T, v, z, dpdv=True)

addon/pyUtils/exports/export_list.py

@@ -188,6 +188,11 @@ def append_export(method, module=None, isBindC=False):
 
 append_export("calc_bmcsl_gij_fmt", "hardsphere_bmcsl")
 
+append_export("hydrate", "init_hydrate_model")
+append_export("hydrate", "fugacity_water_in_hydrate_tpx")
+append_export("hydrate", "fugacity_water_in_hydrate_tvn")
+append_export("hydrate_curves", "map_hydrate_appearance_curve")
+
 append_export("set_standard_entropy", "ideal")
 append_export("get_standard_entropy", "ideal")
 append_export("set_enthalpy_of_formation", "ideal")
@@ -211,6 +216,8 @@ def append_export(method, module=None, isBindC=False):
 append_export("ljs_wca_set_pure_params", "lj_splined")
 append_export("ljs_wca_get_pure_params", "lj_splined")
 
+append_export("multi_phase_envelope_plot_tv", "multi_phase_envelope_tv")
+
 append_export("solve_mu_t", "mut_solver")
 append_export("solve_lnf_t", "mut_solver")
 append_export("map_meta_isotherm", "mut_solver")
@@ -302,6 +309,8 @@ def append_export(method, module=None, isBindC=False):
 append_export("envelope_isentrope_cross", "saturation_curve")
 append_export("pure_fluid_saturation_wrapper", "saturation_curve")
 
+append_export("envelope_plot_tv", "saturation_tv")
+
 append_export("locate_saturation_property", "saturation_point_locators")
 append_export("property_index_from_string", "saturation_point_locators")
 append_export("sat_points_based_on_prop", "saturation_point_locators")

addon/pycThermopack/thermopack/cpa.py

@@ -1,14 +1,7 @@
-# Support for python2
-from __future__ import print_function
-# Import ctypes
 from ctypes import *
-# Importing Numpy (math, arrays, etc...)
 import numpy as np
-# Import platform to detect OS
 from sys import platform, exit
-# Import os utils
 from os import path
-# Import thermo
 from .thermo import c_len_type
 from .cubic import cubic
 

addon/pycThermopack/thermopack/cubic.py

@@ -2,8 +2,8 @@
 import numpy as np
 import warnings
 from os import path
-from .volume_translation import volume_translation
 from .thermo import c_len_type
+from .volume_translation import volume_translation
 
 
 class cubic(volume_translation):