# SWIM Example¶

## Standard Water Integrated Model¶

from thermoengine import phases
from thermoengine import model
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

- H2014, Holten et al., 2014 - W2002, Wagner et al., 2002 - ZD2005, Zhang and Duan, 2005 - DZ2006, Duan and Zhang, 2006
These models are merged, averaged and reconciled for all derivatives to third order over the following regions of $$T$$-$$P$$ space.
plt.plot([198.15,2000.0],[1100.0,1100.0],'y')
plt.plot([198.15,2000.0],[ 900.0, 900.0],'y')
plt.plot([298.15,2000.0],[1000.0,1000.0],'b')
plt.plot([198.15,198.15],[0.0,5000.0],'y')
plt.plot([398.15,398.15],[0.0,5000.0],'y')
plt.plot([298.15,298.15],[0.0,5000.0],'b')
plt.plot([573.15,573.15],[0.0,1100.0],'y')
plt.plot([773.15,773.15],[0.0,1100.0],'y')
plt.plot([673.15,673.15],[0.0,1000.0],'b')
plt.ylabel('P bars')
plt.xlabel('T K')
plt.xlim(left=0.0, right=2000.0)
plt.ylim(bottom=0.0, top=5000.0)
plt.text(1000.0,3000.0,"ZD2005",fontsize=12,bbox=dict(facecolor='orange', alpha=0.2))
plt.text(1200.0, 400.0,"DZ2006",fontsize=12,bbox=dict(facecolor='orange', alpha=0.2))
plt.text( 360.0, 400.0,"W2002",fontsize=12,bbox=dict(facecolor='orange', alpha=0.2))
plt.text(  25.0,2500.0,"H2014",fontsize=12,bbox=dict(facecolor='orange', alpha=0.2))
plt.show()


## Instantiate SWIM using the simple python wrappers¶

(Molecular weight in grams/mole)

modelDB = model.Database()
H2O = modelDB.get_phase('H2O')

print (H2O.props['phase_name'])
print (H2O.props['formula'][0])
print (H2O.props['molwt'][0])

Water
H2O
18.0152


## Use the Python wrapper functions to obtain thermodynamic properties of water¶

$$T$$ (temperature, first argument) is in K, and $$P$$ (pressure, second argument) is in bars.

print ("{0:>10s}{1:15.2f}{2:<20s}".format("G", H2O.gibbs_energy(1000.0, 1000.0), 'J/mol'))
print ("{0:>10s}{1:15.2f}{2:<20s}".format("H", H2O.enthalpy(1000.0, 1000.0), 'J/mol'))
print ("{0:>10s}{1:15.2f}{2:<20s}".format("S", H2O.entropy(1000.0, 1000.0), 'J/K-mol'))
print ("{0:>10s}{1:15.3f}{2:<20s}".format("Cp", H2O.heat_capacity(1000.0, 1000.0), 'J/K-mol'))
print ("{0:>10s}{1:15.6e}{2:<20s}".format("dCp/dT", H2O.heat_capacity(1000.0, 1000.0, deriv={'dT':1}), 'J/-K^2-mol'))
print ("{0:>10s}{1:15.3f}{2:<20s}".format("V", H2O.volume(1000.0, 1000.0, deriv={'dT':1}), 'J/bar-mol'))
print ("{0:>10s}{1:15.6e}{2:<20s}".format("dV/dT", H2O.volume(1000.0, 1000.0, deriv={'dT':1}), 'J/bar-K-mol'))
print ("{0:>10s}{1:15.6e}{2:<20s}".format("dv/dP", H2O.volume(1000.0, 1000.0, deriv={'dP':1}), 'J/bar^2-mol'))
print ("{0:>10s}{1:15.6e}{2:<20s}".format("d2V/dT2", H2O.volume(1000.0, 1000.0, deriv={'dT':2}), 'J/bar-K^2-mol'))
print ("{0:>10s}{1:15.6e}{2:<20s}".format("d2V/dTdP", H2O.volume(1000.0, 1000.0, deriv={'dT':1, 'dP':1}), 'J/bar^2-K-mol'))
print ("{0:>10s}{1:15.6e}{2:<20s}".format("d2V/dP2", H2O.volume(1000.0, 1000.0, deriv={'dP':2}), 'J/bar^3-mol'))

       G     -323482.97J/mol
H     -225835.42J/mol
S         167.15J/K-mol
Cp         72.914J/K-mol
dCp/dT  -1.242746e-01J/-K^2-mol
V          0.014J/bar-mol
dV/dT   1.431826e-02J/bar-K-mol
dv/dP  -7.243395e-03J/bar^2-mol
d2V/dT2  -1.818046e-05J/bar-K^2-mol
d2V/dTdP  -1.299463e-05J/bar^2-K-mol
d2V/dP2   1.744388e-05J/bar^3-mol


## Plot the density of water as a function of $$T$$ …¶

… for 20 isobars at 100 to 2000 bars

P_array = np.linspace(100.0, 2000.0, 20, endpoint=True) # 100->2000, 10 bars
T_array = np.linspace(250.0, 1200.0, 100, endpoint=True) # 250->1200,100 K
MW = H2O.props['molwt']
for P in P_array:
Den_array = MW/H2O.volume(T_array, P)/10.0 ## cc
if P < 1000.0:
plt.plot(T_array, Den_array, 'r-', label=str(int(P)))
else:
plt.plot(T_array, Den_array, 'g-', label=str(int(P))+"-ZD")
plt.plot([673.15,673.15],[0.0,1.1],'b')
plt.plot([298.15,298.15],[0.0,1.1],'b')
plt.ylabel('density g/cc')
plt.xlabel('T K')
plt.title("H2014 -> W2002 -> DZ2006")
plt.legend()
fig = plt.gcf()
fig.set_size_inches(11,8)
plt.show()


## Plot the density of water as a function of $$P$$ …¶

… for 11 isotherms at 0 to 100 °C

T_array = np.linspace(0.0, 100.0, 11, endpoint=True) # 0->100, 11 °C
P_array = np.linspace(100.0, 1000.0, 100, endpoint=True) # 100->1000, 100 bars
MW = H2O.props['molwt']
for T in T_array:
Den_array = MW/H2O.volume(T+273.15, P_array)/10.0 ## cc
if T <= 25.0:
plt.plot(P_array, Den_array, 'b-', label=str(int(T)))
elif T <= 400.0:
plt.plot(P_array, Den_array, 'r-', label=str(int(T)))
else:
plt.plot(P_array, Den_array, 'g-', label=str(int(T)))
plt.plot([1000.0,1000.0],[0.9,1.3])
plt.ylabel('density g/cc')
plt.xlabel('P bars')
plt.ylim(bottom=0.9, top=1.1)
plt.title("W2002 (blue) -> W2002 (red) -> DZ2006 (green)")
plt.legend()
fig = plt.gcf()
fig.set_size_inches(11,8)
plt.show()


## Use direct calls to Objective-C functions (via Rubicon) to select specific water models¶

from ctypes import cdll
from ctypes import util
from rubicon.objc import ObjCClass, objc_method
Water = ObjCClass('GenericH2O')
water = Water.alloc().init()

#water.forceModeChoiceTo_("MELTS H2O-CO2 from Duan and Zhang (2006)")
#water.forceModeChoiceTo_("DEW H2O from Zhang and Duan (2005)")
#water.forceModeChoiceTo_("Supercooled H2O from Holten et al. (2014)")
#water.forceModeChoiceTo_("Steam Properties from Wagner et al. (2002)")
water.forceModeChoiceAutomatic()
T_array = np.linspace(0.0, 130.0, 14, endpoint=True)
P_array = np.linspace(100.0, 2000.0, 100, endpoint=True) # bars
MW = H2O.props['molwt']
for T in T_array:
Den_array = np.empty_like(P_array)
i = 0
for P in P_array:
Den_array[i] = MW/water.getVolumeFromT_andP_(T+273.15, P)/10.0 ## cc
i = i + 1
if T <= 25.0:
plt.plot(P_array, Den_array, 'b-', label=str(int(T)))
elif T <= 400.0:
plt.plot(P_array, Den_array, 'r-', label=str(int(T)))
else:
plt.plot(P_array, Den_array, 'g-', label=str(int(T)))
fig = plt.gcf()
fig.set_size_inches(11,8)
plt.plot([1000.0,1000.0],[0.5,1.3])
plt.ylabel('density g/cc')
plt.xlabel('P bars')
plt.ylim(bottom=0.9, top=1.1)
plt.legend()
plt.show()


# Calculations on the steam saturation curve¶

import pandas as pd
def vol(T=25, P=1):
return H2O.volume(T+273.15, P)*10

from scipy.optimize import curve_fit
def func(x, a, b, c, d, e):
return a + b*x + c*x*x + d*x*x*x + e*x*x*x*x

popt, pcov = curve_fit(func, psat_df['T'], psat_df['P'])

popt

array([ 1.44021565e+00, -2.75944904e-02,  3.50602876e-04, -2.44834016e-06,
1.57085668e-08])

fig = plt.gcf()
fig.set_size_inches(11,8)
plt.title('PSAT curve')
plt.plot(psat_df[['T']], psat_df[['P']], "b-")
plt.plot(psat_df[['T']], func(psat_df['T'], *popt), "r-")
plt.ylabel('Pressure, bars')
plt.xlabel('Temperature, $^\\circ$C')
plt.show()

# create Psat line
volume_psat = vol(psat_df[['T']], psat_df[['P']]+1) # increase psat pressure by 1 bar to ensure liquid H2O
plt.plot(psat_df[['T']], volume_psat, "b-", label="SWIM")
Vol_array = np.empty_like(psat_df['T'])

i = 0
water.forceModeChoiceTo_("MELTS H2O-CO2 from Duan and Zhang (2006)")
for t,p in zip(psat_df['T'], psat_df['P']):
Vol_array[i] = water.getVolumeFromT_andP_(t+273.15, p+1.0)*10.0
i = i + 1
plt.plot(psat_df[['T']], Vol_array, "g-", label="DZ2006")

i = 0
water.forceModeChoiceTo_("Steam Properties from Wagner et al. (2002)")
for t,p in zip(psat_df['T'], psat_df['P']):
Vol_array[i] = water.getVolumeFromT_andP_(t+273.15, p+1.0)*10.0
i = i + 1
plt.plot(psat_df[['T']], Vol_array, "m-", label="Wagner")

def func(x, a, b, c, d, e, f):
return a + b*x + c*x*x + d*x*x*x + e*x*x*x*x + f/(x-374.0)
popt, pcov = curve_fit(func, psat_df['T'], Vol_array)
print (popt)
plt.plot(psat_df[['T']], func(psat_df['T'], *popt), "y--")

water.forceModeChoiceAutomatic()

# create 500 bar line
temps = np.arange(0, 1010, 10)
plt.plot(temps, vol(T=temps, P=500), "r-", label="500 bars")

# plot options
fig = plt.gcf()
fig.set_size_inches(11,8)
plt.title('Partial molal volume of water at PSAT+1 bar and 500 bars')
plt.ylabel('V, $cm^{3}\\cdot mol^{-1}$')
plt.xlabel('Temperature, $^\\circ$C')
plt.margins(x=0) # no margins on x axis
plt.ylim([15, 50])
plt.xlim([0, 380])
plt.xticks(np.arange(0, 380, step=50))
plt.legend()
plt.show()

[ 1.84252342e+01 -3.06586710e-02  5.65750627e-04 -2.69937313e-06
4.67555414e-09 -8.89632469e+00]