Ratio between concentrations in different materials

Hello! Could you please tell me how the concentration jump between different materials is calculated in festim? I expected something like
Ca = Sa/Sb * Cb,
where a, b - different materials, C - concentration at the material boundary, S - “a” and “b” solubility.
Or another ratio is used? I specify only D_0 and E_D values for 2 different materials and festim still solves the problem (without S_0 and E_S).

Hi @cruiseraurora

I specify only D_0 and E_D values for 2 different materials and festim still solves the problem (without S_0 and E_S).

Do you set chemical_pot=True in your model settings to account for the chemical potential conservation at interfaces? When this flag is True, you will have to provide either Sieverts’ or Henry’s constant for each material.

Could you please tell me how the concentration jump between different materials is calculated in festim?

Please, check the corresponding theory section for additional information or the paper about the implementation by @remidm et al.

Please, let me know if it helps

I try to set
chemical_pot = True
and I have the same result (and set S_0 and E_S for each material). But why does it work even without that flag?

Could you, please, provide your input script for further assistance?

Running FESTIM v1.3.1 without solubilities and with chemical_pot=True produces errors.

  File "/home/vvkulagin/anaconda3/envs/festim-env/lib/python3.11/site-packages/festim/materials/materials.py", line 325, in solubility_as_function
    F += mat.S_0 * f.exp(-mat.E_S / k_B / T) * vS * dx(mat.id)
                         ^^^^^^^^
TypeError: bad operand type for unary -: 'NoneType'

Consider this MWE:

import festim as F
import numpy as np

model = F.Simulation()

model.mesh = F.MeshFromVertices(vertices=np.linspace(0,1,num=100))

# Define material properties
mat1 = F.Material(
    id=1,
    D_0=1,
    E_D=1,
    #E_S=1,
    #S_0=1,
    borders=[0,0.5]
)
mat2 = F.Material(
    id=2,
    D_0=2,
    E_D=2,
    #E_S=2,
    #S_0=2,
    borders=[0.5,1.0]
)

model.materials = F.Materials([mat1, mat2])
model.T = 300

model.settings = F.Settings(
    absolute_tolerance=1e12,
    relative_tolerance=1e-8,
    transient=False,
    chemical_pot=True
)

model.initialise()
model.run()

You can uncomment solubilities to run the script without errors.

Sorry, I meant that I get the reliable calculation results (TDS spectrum) without solubilities and without a “chemical_pot=True” in my code. In this case I don’t understand how festim calculates concentration jumps at the boundaries of materials without knowing their solubilities

If you don’t set chemical_pot=True then FESTIM will assume continuity of concentration at interfaces and the solubilities in Material will simply be ignored.

The reason why we have chemical_pot defaulting to False is because it is - slightly - less efficient and not needed for one-material simulations.

1 Like

small remark
@cruiseraurora if you observe “jumps” in the concentration field, I guess, it is just the result of solving the diffusion equation with different diffusivities in materials

I just wanted to understand how concentration calculations work in different materials in festim. Now it is clear to me, thank you both!

2 Likes