Mass transport BC assumes continuity of concentration

[RemDelaporteMathurin]

As discussed with @SamueleMeschini the current mass transport boundary condition assumes continuity of concentration at the fluid-solid interface.

The class MassFlux should therefore be modified to account for this discontinuity.

Instead of having:

$$-D \ \nabla(c) \cdot n = h_{mass} \ (c_\infty - c_\mathrm{surface})$$

the formula should be:

$$-D \ \nabla(c) \cdot n = h_{mass} \ (c_\infty - c_\mathrm{fluid, interface})$$

where $c_\mathrm{fluid, interface}/K_H = (c_\mathrm{surface}/K_S)^2$ in the case of Henry-Sievert interface or $c_\mathrm{fluid, interface}/K_{S,1} = c_\mathrm{surface}/K_{S,2}$ in the case of a Sievert-Sievert interface.

The class should account for both cases. Therefore users should provide:

• solubility law of both materials
• Henry/Sievert constant of both materials

[RemDelaporteMathurin]

I think something like this would work. Based on something @MichaelChristenson proposed

For a Sieverts-Sieverts or Henry-Henry interface:

import fenics as f

class SievertsSievertsMassFlux(F.FluxBC):
    def __init__(self, surfaces, h_coeff, c_ext, S_0, E_S, S_ext_0, E_S_ext) -> None:
        self.h_coeff = h_coeff
        self.c_ext = c_ext
        self.S_0 = S_0
        self.E_S= E_S
        self.S_ext_0 = S_ext_0 
        self.E_S_ext = E_S_ext 
        super().__init__(surfaces=surfaces, field="solute")

    def create_form(self, T, solute):
        h_coeff_expr = Expression(sp.printing.ccode(self.h_coeff), t=0, degree=1)
        c_ext_expr = Expression(sp.printing.ccode(self.c_ext), t=0, degree=1)

        S = self.S_0 * f.exp(-self.E_S/ F.k_B / T)
        S_ext = self.S_ext_0 * f.exp(-self.E_S_ext / F.k_B / T)

        c_interface = S_ext * solute / S

        self.form = -h_coeff_expr * (c_interface - c_ext_expr)
        self.sub_expressions = [h_coeff_expr, c_ext_expr]

For a Sieverts-Henry(ext) interface:

(c/K_S)^2 = c/K_H

import fenics as f

class SievertsHenryMassFlux(F.FluxBC):
    def __init__(self, surfaces, h_coeff, c_ext, S_0, E_S, S_ext_0, E_S_ext) -> None:
        self.h_coeff = h_coeff
        self.c_ext = c_ext
        self.S_0 = S_0
        self.E_S= E_S
        self.S_ext_0 = S_ext_0 
        self.E_S_ext = E_S_ext 
        super().__init__(surfaces=surfaces, field="solute")

    def create_form(self, T, solute):
        h_coeff_expr = Expression(sp.printing.ccode(self.h_coeff), t=0, degree=1)
        c_ext_expr = Expression(sp.printing.ccode(self.c_ext), t=0, degree=1)

        S = self.S_0 * f.exp(-self.E_S/ F.k_B / T)
        S_ext = self.S_ext_0 * f.exp(-self.E_S_ext / F.k_B / T)

        c_interface = S_ext  * (solute/ S)**2

        self.form = -h_coeff_expr * (c_interface - c_ext_expr)
        self.sub_expressions = [h_coeff_expr, c_ext_expr]

We could merge these two classes by doing something like

class MassFlux(F.FluxBC):
    def __init__(self, surfaces, h_coeff, c_ext, S_0, E_S, S_ext_0, E_S_ext, law_domain, law_ext) -> None:
        ...

    def create_form(self, T, solute):
        ....
        if self.law_domain == self.law_ext:
            c_interface = S_ext  * solute/ S
        elif self.law_domain == "sieverts":
            c_interface = S_ext  * (solute/ S)**2
        elif self.law_domain == "henry":
            c_interface = S_ext  * (solute/ S)**0.5  # <--- I don't like this power 0.5 here