Source code for mtnlion.models.thermal

"""
Isothermal model extended with thermal modeling
"""

import dolfin as fem
from ufl.core.expr import Expr  # type: ignore

from mtnlion.formula import Arguments, Formula, Variable
from mtnlion.formulas.approximation import LagrangeMultiplier
from mtnlion.models.isothermal import Isothermal

# pylint: disable=invalid-name


class Thermal(Isothermal):
    """
    The thermal model extends the isothermal model to allow the modeling of internal heat generation which is used to
    determine the temperature of the cell at any location.
    """

    def __init__(self, Ns):
        self.Ns = Ns

        super(Thermal, self).__init__(Ns)

        self.rename_parameter("Ds_ref", "Ds")

        self.variables += [
            Variable("T", ["anode", "separator", "cathode"]),
            Variable("lm_T", ["anode-separator", "separator-cathode"]),
        ]

        self.formulas += [
            self.Temperature(),
            self.HeatGeneration(),
            self.HeatGenerationChemical(),
            self.HeatGenerationEntropy(),
            self.JouleHeatingSolid(),
            self.JouleHeatingElectrolyte1(),
            self.JouleHeatingElectrolyte2(),
            self.Ds(),
            self.DeEff(),
            self.SigmaEff(),
            LagrangeMultiplier(["anode-separator", "separator-cathode"], "lm_T", "T"),
        ]

    class AdaptT(Formula):
        """
        An adapter formula to allow existing formulas to use the temperature variable T as if it were still a parameter.
        """

        def __init__(self):
            super(Thermal.AdaptT, self).__init__(name="T", domains=["anode", "separator", "cathode"])
            self.Variables = self.typedef("Variables", "T")

        def form(self, arguments: Arguments, domain: str) -> Expr:
            (T,) = arguments.variables

            return T

    class ExchangeCurrentDensity(Isothermal.ExchangeCurrentDensity):
        """
        The exchange current density is the current in the absence of net electrolysis and at zero overpotential.
        """

        def __init__(self):
            super(Thermal.ExchangeCurrentDensity, self).__init__()
            self.Variables = self.append_arguments(self.Variables, ("T",))
            self.Parameters = self.append_arguments(self.Parameters, "Eact_k, R, Tref".split(", "))

        def form(self, arguments, domain):
            (T,) = arguments.pop_variables(["T"])
            Eact_k, R, Tref = arguments.pop_parameters("Eact_k, R, Tref".split(", "))

            i0 = super(Thermal.ExchangeCurrentDensity, self).form(arguments, domain)

            return i0 * fem.exp(Eact_k / R * (1 / Tref - 1 / T))

    class Temperature(Formula):
        """
        Temperature of the cell.
        """

        def __init__(self):
            super(Thermal.Temperature, self).__init__(domains=["anode", "separator", "cathode"])
            self.Variables = self.typedef("Variables", "T")
            self.Parameters = self.typedef("Parameters", "rho, cp, L, lambd, q")
            self.TimeDiscretization = self.typedef("TimeDiscretization", "dt_T")

        def form(self, arguments, domain):
            (T,) = arguments.variables
            rho, cp, L, lambd, q = arguments.parameters
            (dt,) = arguments.time_discretization

            lhs = rho * cp * L * dt(T.trial()) * T.test()
            rhs = -lambd / L * fem.dot(fem.grad(T.trial(0)), fem.grad(T.test())) + q * L * T.test()

            return lhs - rhs

    class HeatGenerationChemical(Formula):
        """
        Irreversible heat generation due to chemical reactions for each chemical reaction at the interface.
        """

        def __init__(self):
            super(Thermal.HeatGenerationChemical, self).__init__(name="qi", domains=["anode", "cathode"])
            self.Variables = self.typedef("Variables", "j")
            self.Parameters = self.typedef("Parameters", "a_s, F, eta")

        def form(self, arguments, domain):
            (j,) = arguments.variables
            a_s, F, eta = arguments.parameters

            return a_s * F * eta * j

    class HeatGenerationEntropy(Formula):
        """
        Reversible heat generation due to a change in entropy for each chemical reaction at the interface
        """

        def __init__(self):
            super(Thermal.HeatGenerationEntropy, self).__init__(name="qr", domains=["anode", "cathode"])
            self.Variables = self.typedef("Variables", "j, T")
            self.Parameters = self.typedef("Parameters", "a_s, F, dUocp")

        def form(self, arguments, domain):
            j, T = arguments.variables
            a_s, F, dUocp = arguments.parameters

            return a_s * F * j * T * dUocp

    class JouleHeatingSolid(Formula):
        """
        Joule heating due to electrical potential gradient in the solid.
        """

        def __init__(self):
            super(Thermal.JouleHeatingSolid, self).__init__(name="qs", domains=["anode", "cathode"])
            self.Variables = self.typedef("Variables", "phi_s")
            self.Parameters = self.typedef("Parameters", "L, sigma_eff")

        def form(self, arguments, domain):
            (phis,) = arguments.variables
            (L, sigma_eff) = arguments.parameters

            return sigma_eff * fem.dot(fem.grad(phis.trial(0)), fem.grad(phis.trial(0))) / L ** 2

    class JouleHeatingElectrolyte1(Formula):
        """
        Joule heating due to electrical potential gradient in the electrolyte
        """

        def __init__(self):
            super(Thermal.JouleHeatingElectrolyte1, self).__init__(
                name="qe1", domains=["anode", "separator", "cathode"]
            )
            self.Variables = self.typedef("Variables", "phi_e")
            self.Parameters = self.typedef("Parameters", "L, kappa_eff")

        def form(self, arguments, domain):
            (phie,) = arguments.variables
            (L, kappa_eff) = arguments.parameters

            return kappa_eff * fem.dot(fem.grad(phie.trial(0)), fem.grad(phie.trial(0))) / L ** 2

    class JouleHeatingElectrolyte2(Formula):
        """
        Joule heating due to electrical potential gradient in the electrolyte
        """

        def __init__(self):
            super(Thermal.JouleHeatingElectrolyte2, self).__init__(
                name="qe2", domains=["anode", "separator", "cathode"]
            )
            self.Variables = self.typedef("Variables", "phi_e, c_e")
            self.Parameters = self.typedef("Parameters", "L, kappa_Deff")

        def form(self, arguments, domain):
            phie, ce = arguments.variables
            (L, kappa_Deff) = arguments.parameters

            return kappa_Deff * fem.dot(fem.grad(fem.ln(ce.trial(0))), fem.grad(phie.trial(0))) / L ** 2

    class HeatGeneration(Formula):
        """
        Total heat generated in the cell.
        """

        def __init__(self):
            super(Thermal.HeatGeneration, self).__init__(name="q", domains=["anode", "separator", "cathode"])
            self.Parameters = self.typedef("Parameters", "qi, qr, qs, qe1, qe2")

        def form(self, arguments, domain):
            qi, qr, qs, qe1, qe2 = arguments.parameters

            if domain == "separator":
                qi = qr = qs = 0

            return qi + qr + qs + qe1 + qe2

    class KappaEff(Isothermal.KappaEff):
        """
        Effective conductivity of the electrolyte.
        """

        def __init__(self):
            super(Thermal.KappaEff, self).__init__()
            self.Variables = self.append_arguments(self.Variables, ("T",))
            self.Parameters = self.append_arguments(self.Parameters, "Eact_kappa, R, Tref".split(", "))

        def form(self, arguments, domain):
            (T,) = arguments.pop_variables(("T",))
            Eact_kappa, R, Tref = arguments.pop_parameters("Eact_kappa, R, Tref".split(", "))

            kappa_eff = super(Thermal.KappaEff, self).form(arguments, domain)

            return kappa_eff * fem.exp(Eact_kappa / R * (1 / Tref - 1 / T))

    class SigmaEff(Formula):
        """
        Effective conductivity (electrode-dependent parameter), represents a volume averaged conductivity of the solid
        matrix in a porous media in the vicinity of a given point.
        """

        def __init__(self):
            super(Thermal.SigmaEff, self).__init__(name="sigma_eff", domains=["anode", "cathode"])

            self.Variables = self.append_arguments(self.Variables, ("T",))
            self.Parameters = self.append_arguments(
                self.Parameters, "sigma_ref, Eact_sigma, eps_s, brug_sigma, R, Tref".split(", ")
            )

        def form(self, arguments, domain):
            (T,) = arguments.variables
            sigma_ref, Eact_sigma, eps_s, brug_sigma, R, Tref = arguments.parameters

            return sigma_ref * eps_s ** brug_sigma * fem.exp(Eact_sigma / R * (1 / Tref - 1 / T))

    class KappaDEff(Formula):
        """
        kappa_d effective.
        """

        def __init__(self):
            super(Thermal.KappaDEff, self).__init__(name="kappa_Deff", domains=["anode", "separator", "cathode"])
            self.Parameters = self.typedef("Parameters", "kappa_D, kappa_eff")

        def form(self, arguments, domain):
            kappa_D, kappa_eff = arguments.parameters

            return kappa_D * kappa_eff

    class Ds(Formula):
        """
        Solid diffusivity.
        """

        def __init__(self):
            super(Thermal.Ds, self).__init__(name="Ds", domains=["anode", "cathode"])
            self.Variables = self.typedef("Variables", "T")
            self.Parameters = self.typedef("Parameters", "Ds_ref, Eact_Ds, R, Tref")

        def form(self, arguments, domain):
            (T,) = arguments.variables
            Ds_ref, Eact_Ds, R, Tref = arguments.parameters

            return Ds_ref * fem.exp(Eact_Ds / R * (1 / Tref - 1 / T))

    class DeEff(Formula):
        """
        Effective diffusivity of the electrolyte.
        """

        def __init__(self):
            super(Thermal.DeEff, self).__init__(name="De_eff", domains=["anode", "separator", "cathode"])

            self.Variables = self.typedef("Variables", "T")
            self.Parameters = self.typedef("Parameters", "De_ref, eps_e, brug_De, Eact_De, R, Tref")

        def form(self, arguments, domain):
            (T,) = arguments.variables
            De_ref, eps_e, brug_De, Eact_De, R, Tref = arguments.parameters

            return De_ref * eps_e ** brug_De * fem.exp(Eact_De / R * (1 / Tref - 1 / T))