Soret_coef as a callable

docs/source/theory.rst

@@ -41,9 +41,9 @@ FESTIM can include the Soret effect :cite:`Pendergrass1976,Longhurst1985` (also
 .. math::
     :label: eq_Soret
 
-    J = -D \nabla c_\mathrm{m} - D\frac{Q^* c_\mathrm{m}}{R_g T^2} \nabla T
+    J = -D \nabla c_\mathrm{m} - D\frac{Q^* c_\mathrm{m}}{k_B T^2} \nabla T
 
-where :math:`Q^*` is the Soret coefficient (also called heat of transport) and :math:`R_g` is the gas constant.
+where :math:`Q^*` is the Soret coefficient (also called heat of transport) and :math:`k_B` is the Boltzmann constant.
 
 Conservation of chemical potential at interfaces
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -112,7 +112,7 @@ According to the latter, the rate :math:`k(T)` of a thermally activated process
 
     k(T) = k_0 \exp \left[-\frac{E_k}{k_B T} \right]
 
-where :math:`k_0` is the pre-exponential factor, :math:`E_k` is the process activation energy, :math:`k_B` is the Boltzmann constant, and :math:`T` is the temperature. 
+where :math:`k_0` is the pre-exponential factor, :math:`E_k` is the process activation energy, and :math:`T` is the temperature. 
 
 Heat transfer
 ^^^^^^^^^^^^^^

docs/source/userguide/materials.rst

@@ -35,7 +35,7 @@ Some other parameters are optional and are only required for some types of simul
 * :code:`thermal_cond`: the thermal conductivity in W/m/K
 * :code:`heat_capacity`: the heat capacity in J/kg/K
 * :code:`rho`: the volumetric density in kg/m3
-* :code:`H`: the heat of transport in J/mol that takes a dictionary {"free_enthalpy": …, "entropy": …} so that H = free_enthalpy + entropy*T. For more information see :ref:`Soret effect`.
+* :code:`Q`: the heat of transport in eV. For more information see :ref:`Soret effect`.
 
 --------------------
 Integration with HTM

festim/concentration/mobile.py

@@ -1,4 +1,4 @@
-from festim import Concentration, FluxBC, k_B, R, RadioactiveDecay
+from festim import Concentration, FluxBC, k_B, RadioactiveDecay
 from fenics import *
 
 
@@ -80,10 +80,12 @@ def create_diffusion_form(
                 if mesh.type == "cartesian":
                     F += dot(D * grad(c_0), grad(self.test_function)) * dx
                     if soret:
-                        Q = material.free_enthalpy * T.T + material.entropy
+                        Q = material.Q
+                        if callable(Q):
+                            Q = Q(T.T)
                         F += (
                             dot(
-                                D * Q * c_0 / (R * T.T**2) * grad(T.T),
+                                D * Q * c_0 / (k_B * T.T**2) * grad(T.T),
                                 grad(self.test_function),
                             )
                             * dx

festim/exports/derived_quantities/derived_quantities.py

@@ -79,7 +79,7 @@ def assign_properties_to_quantities(self, materials):
             quantity.D = materials.D
             quantity.S = materials.S
             quantity.thermal_cond = materials.thermal_cond
-            quantity.H = materials.H
+            quantity.Q = materials.Q
 
     def compute(self, t):
         # TODO need to support for soret flag in surface flux

festim/exports/derived_quantities/derived_quantity.py

@@ -15,7 +15,7 @@ def __init__(self, field) -> None:
         self.D = None
         self.S = None
         self.thermal_cond = None
-        self.H = None
+        self.Q = None
         self.data = []
         self.t = []
 

festim/exports/derived_quantities/surface_flux.py

@@ -1,4 +1,4 @@
-from festim import SurfaceQuantity, R
+from festim import SurfaceQuantity, k_B
 import fenics as f
 
 
@@ -22,8 +22,8 @@ def compute(self, soret=False):
             flux += f.assemble(
                 self.prop
                 * self.function
-                * self.H
-                / (R * self.T**2)
+                * self.Q
+                / (k_B * self.T**2)
                 * f.dot(f.grad(self.T), self.n)
                 * self.ds(self.surface)
             )

festim/materials/material.py

@@ -18,9 +18,8 @@ class Material:
             be a function of T. Defaults to None.
         borders (list, optional): The borders of the 1D subdomain.
             Only needed in 1D with several materials. Defaults to None.
-        H (dict, optional): heat of transport (J/mol).
-            {"free_enthalpy": ..., "entropy": ...} so that
-            H = free_enthalpy + entropy*T. Defaults to None.
+        Q (float or callable, optional): heat of transport (eV). Can
+            be a function of T. Defaults to None.
         solubility_law (str, optional): the material's solubility law.
             Can be "henry" or "sievert". Defaults to "sievert".
         name (str, optional): name of the material. Defaults to None.
@@ -37,7 +36,7 @@ def __init__(
         heat_capacity=None,
         rho=None,
         borders=None,
-        H=None,
+        Q=None,
         solubility_law="sievert",
         name=None,
     ) -> None:
@@ -51,10 +50,7 @@ def __init__(
         self.heat_capacity = heat_capacity
         self.rho = rho
         self.borders = borders
-        self.H = H
-        if H is not None:
-            self.free_enthalpy = H["free_enthalpy"]
-            self.entropy = H["entropy"]
+        self.Q = Q
         if solubility_law not in ["henry", "sievert"]:
             raise ValueError(
                 "Acceptable values for solubility_law are 'henry' and 'sievert'"

festim/materials/materials.py

@@ -21,7 +21,7 @@ def __init__(self, materials=[]):
         self.thermal_cond = None
         self.heat_capacity = None
         self.density = None
-        self.H = None
+        self.Q = None
 
     def check_borders(self, size):
         """Checks that the borders of the materials match
@@ -131,7 +131,7 @@ def check_consistency(self):
             "heat_capacity": [],
             "rho": [],
             "borders": [],
-            "H": [],
+            "Q": [],
         }
 
         for attr, value in attributes.items():
@@ -261,8 +261,8 @@ def create_properties(self, vm, T):
             self.thermal_cond = ThermalProp(self, vm, T, "thermal_cond", degree=2)
             self.heat_capacity = ThermalProp(self, vm, T, "heat_capacity", degree=2)
             self.density = ThermalProp(self, vm, T, "rho", degree=2)
-        if self.materials[0].H is not None:
-            self.H = HCoeff(self, vm, T, degree=2)
+        if self.materials[0].Q is not None:
+            self.Q = ThermalProp(self, vm, T, "Q", degree=2)
 
     def solubility_as_function(self, mesh, T):
         """
@@ -361,21 +361,3 @@ def eval_cell(self, value, x, ufc_cell):
 
     def value_shape(self):
         return ()
-
-
-class HCoeff(f.UserExpression):
-    def __init__(self, materials, vm, T, **kwargs):
-        super().__init__(kwargs)
-        self._T = T
-        self._vm = vm
-        self._materials = materials
-
-    def eval_cell(self, value, x, ufc_cell):
-        cell = f.Cell(self._vm.mesh(), ufc_cell.index)
-        subdomain_id = self._vm[cell]
-        material = self._materials.find_material_from_id(subdomain_id)
-
-        value[0] = material.free_enthalpy + self._T(x) * material.entropy
-
-    def value_shape(self):
-        return ()

test/system/test_system.py

@@ -474,28 +474,14 @@ def test_run_MMS_soret(tmpdir):
     D_0 = 2
     k_B = festim.k_B
     D = D_0 * sp.exp(-E_D / k_B / T)
-    H = -2
-    S = 3
-    R = festim.R
+    Q = lambda T: -2e-5 * T + 3e-5
     f = sp.diff(u, festim.t) - sp.diff(
-        (
-            D
-            * (
-                sp.diff(u, festim.x)
-                + (H * T + S) * u / (R * T**2) * sp.diff(T, festim.x)
-            )
-        ),
+        (D * (sp.diff(u, festim.x) + Q(T) * u / (k_B * T**2) * sp.diff(T, festim.x))),
         festim.x,
     )
 
     def run(h):
-        my_materials = festim.Materials(
-            [
-                festim.Material(
-                    id=1, D_0=D_0, E_D=E_D, H={"free_enthalpy": H, "entropy": S}
-                )
-            ]
-        )
+        my_materials = festim.Materials([festim.Material(id=1, D_0=D_0, E_D=E_D, Q=Q)])
         my_initial_conditions = [
             festim.InitialCondition(field=0, value=u),
         ]

test/unit/test_exports/test_derived_quantities/test_derived_quantities.py

@@ -103,7 +103,7 @@ class TestAssignPropertiesToQuantities:
     my_mats = Materials()
     my_mats.D = f.Function(V)
     my_mats.S = f.Function(V)
-    my_mats.H = f.Function(V)
+    my_mats.Q = f.Function(V)
     my_mats.thermal_cond = f.Function(V)
     T = f.Function(V)
 
@@ -117,9 +117,9 @@ def test_quantities_have_S(self):
         for quantity in self.my_quantities.derived_quantities:
             assert quantity.S == self.my_mats.S
 
-    def test_quantities_have_H(self):
+    def test_quantities_have_Q(self):
         for quantity in self.my_quantities.derived_quantities:
-            assert quantity.H == self.my_mats.H
+            assert quantity.Q == self.my_mats.Q
 
     def test_quantities_have_thermal_cond(self):
         for quantity in self.my_quantities.derived_quantities:

test/unit/test_exports/test_derived_quantities/test_surface_flux.py

@@ -1,4 +1,4 @@
-from festim import SurfaceFlux, R
+from festim import SurfaceFlux, k_B
 import fenics as f
 
 
@@ -32,7 +32,7 @@ class TestCompute:
     ds = f.Measure("ds", domain=mesh, subdomain_data=surface_markers)
     D = f.Constant(2)
     thermal_cond = f.Constant(3)
-    H = f.Constant(4)
+    Q = f.Constant(4)
 
     surface = 1
     n = f.FacetNormal(mesh)
@@ -43,7 +43,7 @@ class TestCompute:
     my_h_flux.n = n
     my_h_flux.ds = ds
     my_h_flux.T = T
-    my_h_flux.H = H
+    my_h_flux.Q = Q
 
     my_heat_flux = SurfaceFlux("T", surface)
     my_heat_flux.D = D
@@ -73,8 +73,8 @@ def test_h_flux_with_soret(self):
         expected_flux += f.assemble(
             self.D
             * self.c
-            * self.H
-            / (R * self.T**2)
+            * self.Q
+            / (k_B * self.T**2)
             * f.dot(f.grad(self.T), self.n)
             * self.ds(self.surface)
         )