Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

added derivatives for the CORK EoS #562

Merged
merged 1 commit into from
Dec 17, 2023
Merged

added derivatives for the CORK EoS #562

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
199 changes: 189 additions & 10 deletions burnman/eos/cork.py
View file
Open in desktop
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
Open in desktop
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
Loading