Feature request: class for MassTransferFlux

[SamueleMeschini]

It would be usefull to have a class to model convective mass transfer.

The structure would be similar to the convective heat transfer flux

The mass flux should be: $flux = h{mass}(c{ext} - c_{in})$

where $h{mass}$ is the mass transfer coefficient [m/s], $c{ext}$ is the external concentration (in the fluid flow) and $c_{in}$ is the internal concentration.

I tried to implement it as a new class but it seems not working. Tested on the example from the workshop.

class MassFlux(F.boundary_conditions.fluxes.flux_bc.FluxBC):

    FluxBC subclass for convective mass transfer
    -D * grad(c) * n = h_mass * (c_ext - c)

    Args:
        h_mass (float or sp.Expr): mass transfer coefficient (m/s)
        c_ext (float or sp.Expr): external concentration (1/m3)
        surfaces (list or int): the surfaces of the BC

    def __init__(self, h_coeff, c_ext, surfaces) -> None:
        self.h_coeff = h_coeff
        self.c_ext = c_ext
        super().__init__(surfaces=surfaces, field="0")

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

        self.form = -h_coeff * (solute - c_ext)
        self.sub_expressions = [h_coeff, c_ext]

And the results for a 1D simulation is the following (left boundary with a MassFlux BC):

which is equal to a noFlux BC applied to the left boundary.

[RemDelaporteMathurin]

Hi @SamueleMeschini thanks for that very interesting suggestion!

Would you have a complete MWE (minimal working example) that I could run to have a look? Like what are the value of h and c_ext you used to produce this image?

I'm pretty sure we can make this work!

[SamueleMeschini]

Sure!

import festim as F
import numpy as np
import matplotlib.pyplot as plt
import fenics
from fenics import plot
import sympy as sp

# Define advective mass transfer bc
class MassFlux(F.boundary_conditions.fluxes.flux_bc.FluxBC):
    """
    FluxBC subclass for advective mass flux
    -D * grad(c) * n = h_mass * (c - c_ext)

    Args:
        h_mass (float or sp.Expr): mass transfer coefficient (m/s)
        c_ext (float or sp.Expr): external concentration
        surfaces (list or int): the surfaces of the BC
    """

    def __init__(self, h_coeff, c_ext, surfaces) -> None:
        self.h_coeff = h_coeff
        self.c_ext = c_ext
        super().__init__(surfaces=surfaces, field="0")

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

        self.form = -h_coeff * (solute - c_ext)
        self.sub_expressions = [h_coeff, c_ext]

my_model = F.Simulation()
F.MeshFromVertices(vertices=[0, 1, 2, 3, 4, 5, 6, 7, 7.5])

my_model.mesh = F.MeshFromVertices(
    vertices=np.linspace(0, 1e-6, num=1001)
)

my_model.materials = F.Material(id=1, D_0=1.9e-7, E_D=0.2)
my_model.T = F.Temperature(500)  # K
my_model.boundary_conditions = [
    F.DirichletBC(
        surfaces=2,
        value=1e15,   # H/m3/s
        field=0
        ),
    MassFlux(
        surfaces=1,
        h_coeff = 1e20,
        c_ext = 0
    )
]

my_model.sources = [F.Source(1e15, volume=1, field=0)]
my_model.settings = F.Settings(
    absolute_tolerance=1e10,
    relative_tolerance=1e-10,
    transient=False 
    )
my_model.initialise()
my_model.run()
hydrogen_concentration = my_model.h_transport_problem.mobile.solution
plot(hydrogen_concentration)
plt.xlabel("x [m]")
plt.ylabel("$C_T$ [atoms/m$^3$]")
print("\n Hydrogen concentration = {} [atoms/m3]".format(my_model.h_transport_problem.mobile.solution.vector().get_local()[-1]))```

[RemDelaporteMathurin]

I'll try to make an even impler case to test things out!

[RemDelaporteMathurin]

@SamueleMeschini

I found the issue(s).

field argument should be 0 and not "0". That made the code not recognise the field. This is a bug and should be fixed (i'll open another issue).

Moreover, the create_form method takes in 2 arguments: T and solute (even if T is not used here).

This code (no source, c=0 at right surface and mass transfer flux on the left surface):

import festim as F
import numpy as np
import matplotlib.pyplot as plt
import fenics
from fenics import plot
import sympy as sp

# Define advective mass transfer bc
class MassFlux(F.FluxBC):
    """
    FluxBC subclass for advective mass flux
    -D * grad(c) * n = h_mass * (c - c_ext)

    Args:
        h_mass (float or sp.Expr): mass transfer coefficient (m/s)
        c_ext (float or sp.Expr): external concentration
        surfaces (list or int): the surfaces of the BC
    """

    def __init__(self, h_coeff, c_ext, surfaces) -> None:
        self.h_coeff = h_coeff
        self.c_ext = c_ext
        super().__init__(surfaces=surfaces, field=0)  # field=0 instead of field="0"

    def create_form(self, T, solute):  # F.FluxBC has T and solute as arguments
        h_coeff = fenics.Expression(sp.printing.ccode(self.h_coeff), t=0, degree=1)
        c_ext = fenics.Expression(sp.printing.ccode(self.c_ext), t=0, degree=1)

        self.form = -h_coeff * (solute - c_ext)
        self.sub_expressions = [h_coeff, c_ext]

my_model = F.Simulation()
F.MeshFromVertices(vertices=[0, 1, 2, 3, 4, 5, 6, 7, 7.5])

my_model.mesh = F.MeshFromVertices(
    vertices=np.linspace(0, 1e-6, num=1001)
)

my_model.materials = F.Material(id=1, D_0=1.9e-7, E_D=0.2)
my_model.T = F.Temperature(500)  # K
my_model.boundary_conditions = [
    F.DirichletBC(
        surfaces=2,
        value=0,   # H/m3/s
        field=0
        ),
    MassFlux(
        surfaces=1,
        h_coeff = 1e20,
        c_ext = 1e15
    )
]

my_model.settings = F.Settings(
    absolute_tolerance=1e10,
    relative_tolerance=1e-10,
    transient=False 
    )
my_model.initialise()

my_model.run()
hydrogen_concentration = my_model.h_transport_problem.mobile.solution

plot(hydrogen_concentration)
plt.xlabel("x [m]")
plt.ylabel("$C_T$ [atoms/m$^3$]")
plt.savefig("concentration.png")
print("\n Hydrogen concentration = {} [atoms/m3]".format(my_model.h_transport_problem.mobile.solution.vector().get_local()[-1]))

Produces the expected concentration field:

[jhdark]

Love this idea for a BC, seems so intuitive and I think it would be very useful for my own work too !

[SamueleMeschini]

Great, thanks @RemDelaporteMathurin! @jhdark if you are interested in the physics behind this kind og BC you can find more on the Incropera, Frank P., et al. Fundamentals of heat and mass transfer. Vol. 6. New York: Wiley, 1996.

[RemDelaporteMathurin]

@SamueleMeschini would you be interested in opening a PR to add this very useful class to FESTIM? I have many application cases in mind for it!

[RemDelaporteMathurin]

Will be added in next release