From 92d23a346354bef56eecfc51a768f99a0766c32b Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Sat, 16 Dec 2023 23:57:38 +0000 Subject: [PATCH] added derivatives for the CORK EoS --- burnman/eos/cork.py | 199 ++++++++++++++++++++++++++++++++-- tests/test_eos_consistency.py | 13 +++ 2 files changed, 202 insertions(+), 10 deletions(-) diff --git a/burnman/eos/cork.py b/burnman/eos/cork.py index c94f98fab..c71954f96 100644 --- a/burnman/eos/cork.py +++ b/burnman/eos/cork.py @@ -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 @@ -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): """ @@ -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): @@ -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): @@ -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): """ @@ -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): """ @@ -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: diff --git a/tests/test_eos_consistency.py b/tests/test_eos_consistency.py index 72ef8f824..1074ade21 100644 --- a/tests/test_eos_consistency.py +++ b/tests/test_eos_consistency.py @@ -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