Skip to content

Commit

Permalink
added derivatives for the CORK EoS
Browse files Browse the repository at this point in the history
  • Loading branch information
@bobmyhill
bobmyhill committed Dec 16, 2023
1 parent 33904c4 commit 92d23a3
Show file tree
Hide file tree
Showing 2 changed files with 202 additions and 10 deletions.
199 changes: 189 additions & 10 deletions burnman/eos/cork.py
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for the Earth and Planetary Sciences
# Copyright (C) 2012 - 2017 by the BurnMan team, released under the GNU
# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit
# for the Earth and Planetary Sciences
# Copyright (C) 2012 - 2023 by the BurnMan team, released under the GNU
# GPL v2 or later.

# TO DO: Correct heat capacity, volume where internal order-disorder is
Expand Down Expand Up @@ -48,7 +49,10 @@ def grueneisen_parameter(self, pressure, temperature, volume, params):
Returns grueneisen parameter [unitless] as a function of pressure,
temperature, and volume.
"""
return 0.0
alpha = self.thermal_expansivity(pressure, temperature, volume, params)
K_T = self.isothermal_bulk_modulus(pressure, temperature, volume, params)
C_V = self.molar_heat_capacity_v(pressure, temperature, volume, params)
return alpha * K_T * volume / C_V

def volume(self, pressure, temperature, params):
"""
Expand Down Expand Up @@ -77,7 +81,27 @@ def isothermal_bulk_modulus(self, pressure, temperature, volume, params):
Returns isothermal bulk modulus [Pa] as a function of pressure [Pa],
temperature [K], and volume [m^3]. EQ 13+2
"""
return 0.0
cork = cork_variables(
params["cork_params"], params["cork_P"], params["cork_T"], temperature
)
dVdP = -constants.gas_constant * temperature / (pressure * pressure) + (
-cork[0]
* cork[1]
/ np.sqrt(temperature)
* (
1.0
/ np.power(
constants.gas_constant * temperature + cork[1] * pressure, 2.0
)
- 4.0
/ np.power(
constants.gas_constant * temperature + 2.0 * cork[1] * pressure, 2.0
)
)
+ 0.5 * cork[2] / np.sqrt(pressure)
+ cork[3]
)
return -volume / dVdP

# calculate the shear modulus as a function of P, V, and T
def shear_modulus(self, pressure, temperature, volume, params):
Expand All @@ -93,14 +117,38 @@ def molar_heat_capacity_v(self, pressure, temperature, volume, params):
"""
Returns heat capacity at constant volume at the pressure, temperature, and volume [J/K/mol].
"""
return 0.0
C_p = self.molar_heat_capacity_p(pressure, temperature, volume, params)
V = self.volume(pressure, temperature, params)
alpha = self.thermal_expansivity(pressure, temperature, volume, params)
K_T = self.isothermal_bulk_modulus(pressure, temperature, volume, params)
return C_p - V * temperature * alpha * alpha * K_T

def thermal_expansivity(self, pressure, temperature, volume, params):
"""
Returns thermal expansivity at the pressure, temperature, and volume [1/K]
Replace -Pth in EQ 13+1 with P-Pth for non-ambient temperature
"""
return 0.0
cork = cork_variables(
params["cork_params"], params["cork_P"], params["cork_T"], temperature
)
dcorkdT = params["dcorkdT"]
RT = constants.gas_constant * temperature
c1P = cork[1] * pressure
dVdT = constants.gas_constant / pressure + (
-dcorkdT[0]
* constants.gas_constant
* np.sqrt(temperature)
/ ((RT + c1P) * (RT + 2.0 * c1P))
+ 0.5
* cork[0]
* constants.gas_constant
* (-2.0 * np.power(c1P, 2.0) + 3.0 * RT * (c1P + RT))
/ np.sqrt(temperature)
/ (np.power(RT + c1P, 2.0) * np.power(RT + 2.0 * c1P, 2.0))
+ dcorkdT[2] * np.sqrt(pressure)
+ dcorkdT[3] * pressure
)
return dVdT / volume

# Heat capacity at ambient pressure
def molar_heat_capacity_p0(self, temperature, params):
Expand All @@ -120,14 +168,64 @@ def molar_heat_capacity_p(self, pressure, temperature, volume, params):
"""
Returns heat capacity at constant pressure at the pressure, temperature, and volume [J/K/mol]
"""
return 0
T_0 = params["T_0"]
P_relative = pressure - params["P_0"]

Cp0 = self.molar_heat_capacity_p0(temperature, params)

if params["cork_T"] == 0:
d2RTlnfdTdT = 0.0
else:
cork = cork_variables(
params["cork_params"], params["cork_P"], params["cork_T"], temperature
)
dcorkdT = params["dcorkdT"]

RT = constants.gas_constant * temperature
c1P = cork[1] * P_relative
logterm = np.log(RT + c1P) - np.log(RT + 2.0 * c1P)
dlogtermdT = constants.gas_constant * (
(1.0 / (RT + c1P) - 1.0 / (RT + 2.0 * c1P))
)
d2logtermdTdT = -(
constants.gas_constant
* constants.gas_constant
* (
(
1.0 / np.power(RT + c1P, 2.0)
- 1.0 / np.power(RT + 2.0 * c1P, 2.0)
)
)
)

a = dcorkdT[0] / (cork[1] * np.sqrt(temperature)) - cork[0] / (
2.0 * cork[1] * temperature * np.sqrt(temperature)
)
dadT = (
-0.5 * dcorkdT[0] / (cork[1] * temperature * np.sqrt(temperature))
- dcorkdT[0] / (2.0 * cork[1] * temperature * np.sqrt(temperature))
+ 3.0
* cork[0]
/ (4.0 * cork[1] * temperature * temperature * np.sqrt(temperature))
)
b = cork[0] / (cork[1] * np.sqrt(temperature))
dbdT = dcorkdT[0] / (cork[1] * np.sqrt(temperature)) - 0.5 * b / temperature
d2RTlnfdTdT = (
dadT * logterm + a * dlogtermdT + dbdT * dlogtermdT + b * d2logtermdTdT
) # Second temperature derivative of Eq. 8 in Holland and Powell, 1991

return Cp0 - temperature * d2RTlnfdTdT

def adiabatic_bulk_modulus(self, pressure, temperature, volume, params):
"""
Returns adiabatic bulk modulus [Pa] as a function of pressure [Pa],
temperature [K], and volume [m^3].
"""
return 0.0
K_T = self.isothermal_bulk_modulus(pressure, temperature, volume, params)
C_p = self.molar_heat_capacity_p(pressure, temperature, volume, params)
C_v = self.molar_heat_capacity_v(pressure, temperature, volume, params)
K_S = K_T * C_p / C_v
return K_S

def gibbs_free_energy(self, pressure, temperature, volume, params):
"""
Expand Down Expand Up @@ -192,12 +290,84 @@ def gibbs_free_energy(self, pressure, temperature, volume, params):
+ RTlnf
)

# calculate P = P(T0) + Pth
def entropy(self, pressure, temperature, volume, params):
"""
Returns the entropy [J/K/mol] as a function of pressure [Pa]
and temperature [K].
"""
T_0 = params["T_0"]
P_relative = pressure - params["P_0"]

intCpoverTdT = (
params["Cp"][0] * np.log(temperature)
+ params["Cp"][1] * temperature
- 0.5 * params["Cp"][2] / np.power(temperature, 2.0)
- 2.0 * params["Cp"][3] / np.sqrt(temperature)
) - (
params["Cp"][0] * np.log(T_0)
+ params["Cp"][1] * T_0
- 0.5 * params["Cp"][2] / (T_0 * T_0)
- 2.0 * params["Cp"][3] / np.sqrt(T_0)
)

if params["cork_T"] == 0:
dRTlnfdT = 0.0
else:
cork = cork_variables(
params["cork_params"], params["cork_P"], params["cork_T"], temperature
)
dcorkdT = params["dcorkdT"]

logterm = np.log(
constants.gas_constant * temperature + cork[1] * P_relative
) - np.log(
constants.gas_constant * temperature + 2.0 * cork[1] * P_relative
)
dRTlnfdT = (
constants.gas_constant * np.log(1e-5 * P_relative)
+ (
dcorkdT[0] / (cork[1] * np.sqrt(temperature))
- cork[0] / (2.0 * cork[1] * temperature * np.sqrt(temperature))
)
* logterm
+ (cork[0] / (cork[1] * np.sqrt(temperature)))
* constants.gas_constant
* (
(
1.0
/ (constants.gas_constant * temperature + cork[1] * P_relative)
- 1.0
/ (
constants.gas_constant * temperature
+ 2.0 * cork[1] * P_relative
)
)
)
+ 2.0 / 3.0 * dcorkdT[2] * P_relative * np.sqrt(P_relative)
+ dcorkdT[3] / 2.0 * P_relative * P_relative
) # Temperature derivative of Eq. 8 in Holland and Powell, 1991

return params["S_0"] + intCpoverTdT - dRTlnfdT

def helmholtz_free_energy(self, pressure, temperature, volume, params):
return self.gibbs_free_energy(
pressure, temperature, volume, params
) - pressure * self.volume(pressure, temperature, params)

def enthalpy(self, pressure, temperature, volume, params):
"""
Returns the enthalpy [J/mol] as a function of pressure [Pa]
and temperature [K].
"""
gibbs = self.gibbs_free_energy(pressure, temperature, volume, params)
entropy = self.entropy(pressure, temperature, volume, params)
return gibbs + temperature * entropy

def pressure(self, temperature, volume, params):
"""
Returns pressure [Pa] as a function of temperature [K] and volume[m^3]
"""
return 0.0
return NotImplementedError("")

def validate_parameters(self, params):
"""
Expand All @@ -217,6 +387,15 @@ def validate_parameters(self, params):
if "S_0" not in params:
params["S_0"] = float("nan")

cork, cork_P, cork_T = params["cork_params"], params["cork_P"], params["cork_T"]
if "dcorkdT" not in params:
params["dcorkdT"] = [
cork[0][1] * cork_T ** (1.5) / cork_P,
0.0,
cork[2][1] / cork_P ** (1.5),
cork[3][1] / cork_P ** (2.0),
]

# check that all the required keys are in the dictionary
expected_keys = ["cork_params", "cork_T", "cork_P", "Cp"]
for k in expected_keys:
Expand Down
13 changes: 13 additions & 0 deletions tests/test_eos_consistency.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,19 @@ def test_HP(self):
True,
)

def test_CORK(self):
P = 10.0e9
T = 3000.0
self.assertEqual(
check_eos_consistency(
burnman.minerals.HP_2011_fluids.CO2(),
P,
T,
including_shear_properties=False,
),
True,
)

def test_SLB(self):
P = 10.0e9
T = 3000.0
Expand Down

0 comments on commit 92d23a3

Please sign in to comment.