Skip to content

Commit c2ee8b6

Browse files
bobmyhillbobmyhill
committed
added simple anharmonic model
1 parent c674004 commit c2ee8b6

File tree

4 files changed

+323
-0
lines changed

4 files changed

+323
-0
lines changed

‎burnman/eos/helper.py‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -63,6 +63,8 @@ def create(method):
6363
return slb.SLB3Conductive()
6464
elif method == "modular_mgd":
6565
return mmgd.ModularMGD()
66+
elif method == "modular_mgd_with_anharmonicity":
67+
return mmgd.ModularMGDWithAnharmonicity()
6668
elif method == "murnaghan":
6769
return murnaghan.Murnaghan()
6870
elif method == "bm3shear2":

‎burnman/eos/modular_mie_grueneisen_debye.py‎

Lines changed: 295 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -20,9 +20,11 @@ def decorator(fn):
2020
return decorator
2121

2222

23+
import numpy as np
2324
from . import debye
2425
from . import equation_of_state as eos
2526
from ..utils.math import bracket
27+
from ..utils.misc import copy_documentation
2628
from . import bukowinski_electronic as el
2729
from .anharmonic_debye import AnharmonicDebye as Anharmonic
2830

@@ -351,3 +353,296 @@ def validate_parameters(self, params):
351353
if "bel_0" not in params:
352354
params["bel_0"] = None
353355
params["gel"] = None
356+
357+
358+
class ModularMGDWithAnharmonicity(ModularMGD):
359+
"""
360+
This class extends the ModularMGD class to include anharmonicity effects
361+
according to a simplification of the model proposed by
362+
Wu and Wentzcovitch, 2009.
363+
364+
The basis of the anharmonic model is the definition of a scaled volume, :math:`V'`:
365+
:math:`\\ln (V'/V) = = - c \\ln (V/V_0)`.
366+
In this expression, :math:`V` is the target volume,
367+
:math:`V_0` is a first order approximation to the volume at the same pressure
368+
and reference temperature :math:`T_0`, and :math:`c` is an anharmonicity parameter
369+
provided in the params dictionary as `c_anh`.
370+
371+
The anharmonic Helmholtz energy :math:`F` is related to the scaled volume by the equation:
372+
373+
.. math::
374+
375+
F(V,T) = F_h(V',T) + F_h(V,T_0) - F_h(V',T_0)
376+
377+
where :math:`F_h` is the harmonic Helmholtz energy,
378+
potentially with electronic contributions.
379+
380+
Note: This model is not the same as that published in
381+
Wu and Wentzcovitch (2009). The results are expected to be similar,
382+
but the :math:`c` parameter will in general need to be tweaked.
383+
This is because only a local approximation to the volume change
384+
between 0 K and the target temperature is used. This does not mean that
385+
the model is less able to capture the essential physics of the problem;
386+
indeed, the model of Wu and Wentzcovitch (2009) is only intended to be
387+
an effective ansatz.
388+
"""
389+
390+
def lnVoverV0_approx(self, temperature, volume, params):
391+
a_h = ModularMGD.thermal_expansivity(
392+
super(), np.nan, temperature, volume, params
393+
)
394+
KT_h = ModularMGD.isothermal_bulk_modulus_reuss(
395+
super(), np.nan, temperature, volume, params
396+
)
397+
KT0_h = ModularMGD.isothermal_bulk_modulus_reuss(
398+
super(), np.nan, params["T_0"], volume, params
399+
)
400+
P_th = a_h * KT_h * (temperature - params["T_0"])
401+
return P_th / KT0_h # assume a constant bulk modulus
402+
403+
@copy_documentation(ModularMGD.pressure)
404+
def pressure(self, temperature, volume, params):
405+
dV = 1e-5 * params["V_0"]
406+
407+
lnVoverV0 = self.lnVoverV0_approx(temperature, volume, params)
408+
dlnVoverV0_dV = (
409+
self.lnVoverV0_approx(temperature, volume + dV, params)
410+
- self.lnVoverV0_approx(temperature, volume - dV, params)
411+
) / (2 * dV)
412+
413+
lnVprimeoverV = -params["c_anh"] * lnVoverV0
414+
Vprime = volume * np.exp(lnVprimeoverV)
415+
416+
dlnVprimeoverV_dV = -params["c_anh"] * dlnVoverV0_dV
417+
dVprime_dV = Vprime * (1.0 / volume + dlnVprimeoverV_dV)
418+
419+
T0 = params["T_0"]
420+
P_V_T0 = ModularMGD.pressure(super(), T0, volume, params)
421+
P_Vp_T = ModularMGD.pressure(super(), temperature, Vprime, params)
422+
P_Vp_T0 = ModularMGD.pressure(super(), T0, Vprime, params)
423+
424+
P = P_V_T0 + dVprime_dV * (P_Vp_T - P_Vp_T0)
425+
426+
return P
427+
428+
@copy_documentation(ModularMGD.isothermal_bulk_modulus_reuss)
429+
def isothermal_bulk_modulus_reuss(self, pressure, temperature, volume, params):
430+
dV = 1e-5 * params["V_0"]
431+
432+
lnVoverV0 = self.lnVoverV0_approx(temperature, volume, params)
433+
lnVoverV0_plus_dV = self.lnVoverV0_approx(temperature, volume + dV, params)
434+
lnVoverV0_minus_dV = self.lnVoverV0_approx(temperature, volume - dV, params)
435+
436+
dlnVoverV0_dV = (lnVoverV0_plus_dV - lnVoverV0_minus_dV) / (2 * dV)
437+
d2lnVoverV0_dV2 = (lnVoverV0_plus_dV - 2 * lnVoverV0 + lnVoverV0_minus_dV) / (
438+
dV**2
439+
)
440+
441+
lnVprimeoverV = -params["c_anh"] * lnVoverV0
442+
Vprime = volume * np.exp(lnVprimeoverV)
443+
444+
dlnVprimeoverV_dV = -params["c_anh"] * dlnVoverV0_dV
445+
dVprime_dV = Vprime * (1.0 / volume + dlnVprimeoverV_dV)
446+
d2Vprime_dV2 = Vprime * (
447+
(1.0 / volume + dlnVprimeoverV_dV) ** 2
448+
- 1.0 / volume**2
449+
- params["c_anh"] * d2lnVoverV0_dV2
450+
)
451+
452+
T0 = params["T_0"]
453+
P_Vp_T = ModularMGD.pressure(super(), temperature, Vprime, params)
454+
P_Vp_T0 = ModularMGD.pressure(super(), T0, Vprime, params)
455+
456+
KT_V_T0 = ModularMGD.isothermal_bulk_modulus_reuss(
457+
super(), np.nan, T0, volume, params
458+
)
459+
KT_Vp_T = ModularMGD.isothermal_bulk_modulus_reuss(
460+
super(), np.nan, temperature, Vprime, params
461+
)
462+
KT_Vp_T0 = ModularMGD.isothermal_bulk_modulus_reuss(
463+
super(), np.nan, T0, Vprime, params
464+
)
465+
466+
KT = (
467+
KT_V_T0
468+
- volume * d2Vprime_dV2 * (P_Vp_T - P_Vp_T0)
469+
+ volume / Vprime * (dVprime_dV) ** 2 * (KT_Vp_T - KT_Vp_T0)
470+
)
471+
return KT
472+
473+
@copy_documentation(ModularMGD._molar_heat_capacity_v)
474+
def _molar_heat_capacity_v(self, pressure, temperature, volume, params):
475+
dT = 1e-2
476+
477+
lnVoverV0 = self.lnVoverV0_approx(temperature, volume, params)
478+
lnVoverV0_plus_dT = self.lnVoverV0_approx(temperature + dT, volume, params)
479+
lnVoverV0_minus_dT = self.lnVoverV0_approx(temperature - dT, volume, params)
480+
481+
dlnVoverV0_dT = (lnVoverV0_plus_dT - lnVoverV0_minus_dT) / (2 * dT)
482+
d2lnVoverV0_dT2 = (lnVoverV0_plus_dT - 2 * lnVoverV0 + lnVoverV0_minus_dT) / (
483+
dT**2
484+
)
485+
486+
c = params["c_anh"]
487+
lnVprimeoverV = -c * lnVoverV0
488+
Vprime = volume * np.exp(lnVprimeoverV)
489+
490+
dlnVprimeoverV_dT = -c * dlnVoverV0_dT
491+
d2lnVprimeoverV_dT2 = -c * d2lnVoverV0_dT2
492+
493+
dVprime_dT = Vprime * dlnVprimeoverV_dT
494+
d2Vprime_dT2 = Vprime * (dlnVprimeoverV_dT**2 + d2lnVprimeoverV_dT2)
495+
496+
T0 = params["T_0"]
497+
498+
Cvh_Vprime_T = ModularMGD._molar_heat_capacity_v(
499+
super(), np.nan, temperature, Vprime, params
500+
)
501+
P_Vprime_T = ModularMGD.pressure(super(), temperature, Vprime, params)
502+
P_Vprime_T0 = ModularMGD.pressure(super(), T0, Vprime, params)
503+
KT_Vprime_T = ModularMGD.isothermal_bulk_modulus_reuss(
504+
super(), np.nan, temperature, Vprime, params
505+
)
506+
KT_Vprime_T0 = ModularMGD.isothermal_bulk_modulus_reuss(
507+
super(), np.nan, T0, Vprime, params
508+
)
509+
dPdT_Vprime_T = (
510+
ModularMGD.thermal_expansivity(super(), np.nan, temperature, Vprime, params)
511+
* KT_Vprime_T
512+
)
513+
514+
C_v = Cvh_Vprime_T + temperature * (
515+
2.0 * dPdT_Vprime_T * dVprime_dT
516+
- (P_Vprime_T0 - P_Vprime_T) * d2Vprime_dT2
517+
- (KT_Vprime_T - KT_Vprime_T0) / Vprime * dVprime_dT**2
518+
)
519+
return C_v
520+
521+
def thermal_expansivity(self, pressure, temperature, volume, params):
522+
dV = 1e-5 * params["V_0"]
523+
dT = 1e-2
524+
525+
# lnVoverV0 derivatives
526+
lnVoverV0 = self.lnVoverV0_approx(temperature, volume, params)
527+
lnVoverV0_plusV = self.lnVoverV0_approx(temperature, volume + dV, params)
528+
lnVoverV0_minusV = self.lnVoverV0_approx(temperature, volume - dV, params)
529+
lnVoverV0_plusT = self.lnVoverV0_approx(temperature + dT, volume, params)
530+
lnVoverV0_minusT = self.lnVoverV0_approx(temperature - dT, volume, params)
531+
lnVoverV0_plusV_plusT = self.lnVoverV0_approx(
532+
temperature + dT, volume + dV, params
533+
)
534+
lnVoverV0_plusV_minusT = self.lnVoverV0_approx(
535+
temperature - dT, volume + dV, params
536+
)
537+
lnVoverV0_minusV_plusT = self.lnVoverV0_approx(
538+
temperature + dT, volume - dV, params
539+
)
540+
lnVoverV0_minusV_minusT = self.lnVoverV0_approx(
541+
temperature - dT, volume - dV, params
542+
)
543+
544+
dlnVoverV0_dV = (lnVoverV0_plusV - lnVoverV0_minusV) / (2 * dV)
545+
dlnVoverV0_dT = (lnVoverV0_plusT - lnVoverV0_minusT) / (2 * dT)
546+
d2lnVoverV0_dTdV = (
547+
lnVoverV0_plusV_plusT
548+
- lnVoverV0_plusV_minusT
549+
- lnVoverV0_minusV_plusT
550+
+ lnVoverV0_minusV_minusT
551+
) / (4 * dV * dT)
552+
553+
c = params["c_anh"]
554+
lnVprimeoverV = -c * lnVoverV0
555+
Vprime = volume * np.exp(lnVprimeoverV)
556+
dlnVprime_dV = -c * dlnVoverV0_dV
557+
dlnVprime_dT = -c * dlnVoverV0_dT
558+
d2lnVprime_dTdV = -c * d2lnVoverV0_dTdV
559+
560+
dVprime_dV = Vprime * (1.0 / volume + dlnVprime_dV)
561+
dVprime_dT = Vprime * dlnVprime_dT
562+
d2Vprime_dTdV = Vprime * (
563+
dlnVprime_dT * (1.0 / volume + dlnVprime_dV) + d2lnVprime_dTdV
564+
)
565+
566+
T0 = params["T_0"]
567+
P_Vprime_T = ModularMGD.pressure(super(), temperature, Vprime, params)
568+
P_Vprime_T0 = ModularMGD.pressure(super(), T0, Vprime, params)
569+
KT_Vprime_T = ModularMGD.isothermal_bulk_modulus_reuss(
570+
super(), np.nan, temperature, Vprime, params
571+
)
572+
KT_Vprime_T0 = ModularMGD.isothermal_bulk_modulus_reuss(
573+
super(), np.nan, T0, Vprime, params
574+
)
575+
alphaKT_Vprime_T = (
576+
ModularMGD.thermal_expansivity(super(), np.nan, temperature, Vprime, params)
577+
* KT_Vprime_T
578+
)
579+
580+
aKT = (
581+
alphaKT_Vprime_T * dVprime_dV
582+
- (KT_Vprime_T - KT_Vprime_T0) * dVprime_dV * dVprime_dT / Vprime
583+
+ (P_Vprime_T - P_Vprime_T0) * d2Vprime_dTdV
584+
)
585+
KT = self.isothermal_bulk_modulus_reuss(pressure, temperature, volume, params)
586+
587+
return aKT / KT
588+
589+
def entropy(self, pressure, temperature, volume, params):
590+
dT = 1e-2
591+
592+
lnVoverV0 = self.lnVoverV0_approx(temperature, volume, params)
593+
dlnVoverV0_dT = (
594+
self.lnVoverV0_approx(temperature + dT, volume, params)
595+
- self.lnVoverV0_approx(temperature - dT, volume, params)
596+
) / (2 * dT)
597+
598+
c = params["c_anh"]
599+
lnVprimeoverV = -c * lnVoverV0
600+
Vprime = volume * np.exp(lnVprimeoverV)
601+
602+
dlnVprimeoverV_dT = -c * dlnVoverV0_dT
603+
dVprime_dT = Vprime * dlnVprimeoverV_dT
604+
605+
Sh_Vprime_T = ModularMGD.entropy(super(), np.nan, temperature, Vprime, params)
606+
P_Vprime_T = ModularMGD.pressure(super(), temperature, Vprime, params)
607+
P_Vprime_T0 = ModularMGD.pressure(super(), params["T_0"], Vprime, params)
608+
609+
S = Sh_Vprime_T + (P_Vprime_T - P_Vprime_T0) * dVprime_dT
610+
return S
611+
612+
def _helmholtz_energy(self, pressure, temperature, volume, params):
613+
lnVoverV0 = self.lnVoverV0_approx(temperature, volume, params)
614+
lnVprimeoverV = -params["c_anh"] * lnVoverV0
615+
Vprime = volume * np.exp(lnVprimeoverV)
616+
Fh_Vprime_T = ModularMGD._helmholtz_energy(
617+
super(), np.nan, temperature, Vprime, params
618+
)
619+
F_V_T0 = ModularMGD._helmholtz_energy(
620+
super(), np.nan, params["T_0"], volume, params
621+
)
622+
F_Vprime_T0 = ModularMGD._helmholtz_energy(
623+
super(), np.nan, params["T_0"], Vprime, params
624+
)
625+
626+
F = Fh_Vprime_T + F_V_T0 - F_Vprime_T0
627+
return F
628+
629+
def validate_parameters(self, params):
630+
"""
631+
Check for existence and validity of the parameters
632+
"""
633+
634+
ModularMGD.validate_parameters(self, params)
635+
636+
if "anharmonic_thermal_model" in params:
637+
raise ValueError(
638+
"The ModularMGDWithAnharmonicity class "
639+
"does not allow for a separate anharmonic "
640+
"contribution. Anharmonicity is incorporated "
641+
"through a single 'c_anh' parameter."
642+
)
643+
644+
if "c_anh" not in params:
645+
raise ValueError(
646+
"The ModularMGDWithAnharmonicity class "
647+
"requires a 'c_anh' parameter to be set."
648+
)

‎docs/eos.rst‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -104,6 +104,8 @@ Modular Mie-Grüneisen-Debye
104104

105105
.. autoclass:: burnman.eos.ModularMGD
106106

107+
.. autoclass:: burnman.eos.ModularMGDWithAnharmonicity
108+
107109
Debye Temperature Models
108110
^^^^^^^^^^^^^^^^^^^^^^^^
109111

‎tests/test_modular_mgd.py‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -281,6 +281,30 @@ def test_SLB_lognormal_anharmonic_contribution(self):
281281
m_without_component.set_state(2.0e9, 2000.0)
282282
self.assertFalse(m_without_component.S == m.S)
283283

284+
def test_SLB_inbuilt_anharmonic_contribution(self):
285+
params = mineral_params.copy()
286+
params["reference_eos"] = create("bm3")
287+
params["debye_temperature_model"] = debye_temperature_models.SLB()
288+
params["Debye_0"] = 1000.0
289+
params["grueneisen_0"] = 1.2
290+
params["q_0"] = 1.1
291+
params["P_0"] = 1.0e5
292+
m_without_component = Mineral(params)
293+
294+
params2 = params.copy()
295+
params2["equation_of_state"] = "modular_mgd_with_anharmonicity"
296+
params2["c_anh"] = 0.05
297+
298+
m = Mineral(params2)
299+
consistent = check_eos_consistency(
300+
m, 2.0e9, 2000.0, including_shear_properties=False, tol=1.0e-4
301+
)
302+
self.assertTrue(consistent)
303+
304+
m.set_state(2.0e9, 2000.0)
305+
m_without_component.set_state(2.0e9, 2000.0)
306+
self.assertFalse(m_without_component.S == m.S)
307+
284308
def test_SPOCK_isothermal_contribution(self):
285309
params = mineral_params.copy()
286310
params["reference_eos"] = create("spock")

0 commit comments

Comments
 (0)