Multi-material hydrogen permeation 3D model transient solution not converging

[jaronfcota]

Hello,

I am having a bit of trouble with a model converging. I’m trying to model hydrogen permeation through a nickel membrane through molten FLiBe in a cylindrical geometry, and when I try and calculate the transient solution to the flux of the top surface of the FLiBe the solver is not converging.

Here is an example picture of the model:

There is a nickel membrane, a layer of FLiBe salt, and a nickel wall surrounding both the membrane and the salt.

Here is an example of the log when running the model:

Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 3.143e+25 (tol = 1.000e+20) r (rel) = 1.000e+00 (tol = 1.000e-10)
  Newton iteration 1: r (abs) = 4.078e+24 (tol = 1.000e+20) r (rel) = 1.298e-01 (tol = 1.000e-10)
  Newton iteration 2: r (abs) = 6.666e+23 (tol = 1.000e+20) r (rel) = 2.121e-02 (tol = 1.000e-10)
  Newton iteration 3: r (abs) = 1.949e+23 (tol = 1.000e+20) r (rel) = 6.200e-03 (tol = 1.000e-10)
  Newton iteration 4: r (abs) = 1.178e+23 (tol = 1.000e+20) r (rel) = 3.747e-03 (tol = 1.000e-10)
  Newton iteration 5: r (abs) = 7.058e+28 (tol = 1.000e+20) r (rel) = 2.246e+03 (tol = 1.000e-10)
  Newton iteration 6: r (abs) = 1.765e+28 (tol = 1.000e+20) r (rel) = 5.614e+02 (tol = 1.000e-10)
  Newton iteration 7: r (abs) = 4.411e+27 (tol = 1.000e+20) r (rel) = 1.404e+02 (tol = 1.000e-10)
  Newton iteration 8: r (abs) = 1.113e+27 (tol = 1.000e+20) r (rel) = 3.540e+01 (tol = 1.000e-10)
  Newton iteration 9: r (abs) = 2.782e+26 (tol = 1.000e+20) r (rel) = 8.852e+00 (tol = 1.000e-10)
  Newton iteration 10: r (abs) = 3.863e+26 (tol = 1.000e+20) r (rel) = 1.229e+01 (tol = 1.000e-10)
  Newton iteration 11: r (abs) = 9.663e+25 (tol = 1.000e+20) r (rel) = 3.075e+00 (tol = 1.000e-10)
  Newton iteration 12: r (abs) = 2.638e+25 (tol = 1.000e+20) r (rel) = 8.393e-01 (tol = 1.000e-10)
  Newton iteration 13: r (abs) = 1.408e+25 (tol = 1.000e+20) r (rel) = 4.480e-01 (tol = 1.000e-10)

Here are the meshing files generated of the geometry

8mm_thick_mesh.zip (910.8 KB)

And these meshing files gives the ids of

{-6: ['salt_volume'], -7: ['salt_top'], -8: ['membrane_bottom'], -9: ['side_wall'], -10: ['nickel_volume']}

And here is the code that I am using to try and solve the transient problem:

import festim as F
import fenics as f
import h_transport_materials as htm
import numpy as np

P_up = 100

def make_model(folder):
    model = F.Simulation()

    model.mesh = F.MeshFromXDMF(
        volume_file=f"{folder}/mesh_domains.xdmf", boundary_file=f"{folder}/mesh_boundaries.xdmf"
    )

    # Material properties
    flibe_diffusivity = htm.diffusivities.filter(material=htm.FLIBE, author = 'calderoni')[0]
    flibe_solubility = htm.solubilities.filter(material=htm.FLIBE, author = 'calderoni')[0]
    flibe = F.Material(
        id=6, # Given by the meshing file
        D_0=flibe_diffusivity.pre_exp.magnitude,
        E_D=flibe_diffusivity.act_energy.magnitude,
        S_0 = flibe_solubility.pre_exp.magnitude,
        E_S = flibe_solubility.act_energy.magnitude,
        solubility_law = "henry"
    )

    nickel_diffusivity = htm.diffusivities.filter(material=htm.NICKEL).mean()
    nickel_solubility = htm.solubilities.filter(material=htm.NICKEL).mean()
    nickel = F.Material(
            id = 10,
            D_0 = nickel_diffusivity.pre_exp.magnitude,
            E_D = nickel_diffusivity.act_energy.magnitude,
            S_0 = nickel_solubility.pre_exp.magnitude,
            E_S = nickel_solubility.act_energy.magnitude
    )

    model.materials = [nickel, flibe]

    model.boundary_conditions = [
        F.SievertsBC(surfaces=[8], S_0=nickel.S_0, E_S=nickel.E_S, pressure = P_up),
        # Dirichlet BCs on salt top and outside of side wall
        F.DirichletBC(field="solute", value=0, surfaces=[7]),
        F.DirichletBC(field="solute", value = 0, surfaces = [9])
    ]

    model.T = F.Temperature(973)

    model.exports = [F.XDMFExport("solute", folder = folder)]

    model.settings = F.Settings(
        absolute_tolerance=1e20,
        relative_tolerance=1e-10,
        chemical_pot=True,
        final_time = 8*3600
    )

    model.dt = F.Stepsize(initial_value = 10, stepsize_change_ratio = 1.1)

    model.log_level = 20

    return model

if __name__ == "__main__":


    folder = '3D_cad/8mm_thick'
    model = make_model(folder)

    model.initialise()
    model.run()

I suspect that the meshing might have something to do with the lack of convergence or maybe an improper implementation of boundary conditions. Some help and guidance would be much appreciated!

[remidm]

Hi @jaronfcota

I haven’t run your model yet but this could be a coarse mesh issue.

Have you tried running the simulation in steady state?

Have you considered running a 2D axisymmetric model instead of a 3D model?

I would simply mesh half of the domain in the (r,z) plane and set the mesh to cylindrical with

MeshFromXDMF(……., type=“cylindrical”)

That way, you can have a much finer mesh and probably better converge to a solution

[remidm]

Looking at the mesh now, I think it may be because the mesh is too coarse. there is only four cells in the z direction.

Maybe this mesh is enough if you solve your problem in steady state by:

• removing the stepsize
• setting the transient argument of F.Settings to False

To have a transient model, I strongly recommend doing it in 2D axisymmetric. Consider this example for non-cartesian meshes.
In SALOME, mesh only a (r,z) cross section (only half the geometry) and you can increase the number of cells.

[jaronfcota]

Yeah that makes sense. This model does solve the steady state problem:

mobile_concentration.xdmf (1.6 KB)

I did notice that there was an example for axisymmetric geometries and that is probably just a better approach overall. I’ll tinker around with refining the mesh and report if that gives transient results. Thanks!

[jaronfcota]

Hi,

After creating a 2D axisymmetric model I’m still getting questionable estimations of the surface flux of the salt.

Again here are my meshing files and SALOME CAD files:
CAD_and_meshes.zip (1.3 MB)

Resulting in a correspondence dictionary of

{-6: ['salt_face'], -7: ['nickel_faces'], -8: ['membrane_bottom'], -9: ['side_wall'], -10: ['salt_top']}

And again here is the code I am using to evaluate this mesh:

import festim as F
import fenics as f
import h_transport_materials as htm
import numpy as np
import matplotlib.pyplot as plt


P_up = 100

def make_model(folder):
    model = F.Simulation()

    model.mesh = F.MeshFromXDMF(
        volume_file=f"{folder}/mesh_domains.xdmf", boundary_file=f"{folder}/mesh_boundaries.xdmf", type = 'cylindrical'
    )

    # Material properties
    # Using calderoni's numbers instead of the mean
    flibe_diffusivity = htm.diffusivities.filter(material=htm.FLIBE, author = 'calderoni')[0]
    # Adding flibe solubility properties
    flibe_solubility = htm.solubilities.filter(material=htm.FLIBE, author = 'calderoni')[0]
    flibe = F.Material(
        id=6, # Given by the meshing file
        D_0=flibe_diffusivity.pre_exp.magnitude,
        E_D=flibe_diffusivity.act_energy.magnitude,
        S_0 = flibe_solubility.pre_exp.magnitude,
        E_S = flibe_solubility.act_energy.magnitude,
        solubility_law = "henry"
    )

    nickel_diffusivity = htm.diffusivities.filter(material=htm.NICKEL).mean()
    nickel_solubility = htm.solubilities.filter(material=htm.NICKEL).mean()
    nickel = F.Material(
            id = 7,
            D_0 = nickel_diffusivity.pre_exp.magnitude,
            E_D = nickel_diffusivity.act_energy.magnitude,
            S_0 = nickel_solubility.pre_exp.magnitude,
            E_S = nickel_solubility.act_energy.magnitude
    )

    model.materials = [nickel, flibe]

    model.boundary_conditions = [
        F.SievertsBC(surfaces=[8], S_0=nickel.S_0, E_S=nickel.E_S, pressure = P_up),
        # Dirichlet BCs on salt top and outside of side wall
        F.DirichletBC(field="solute", value = 0, surfaces=[9]),
        F.DirichletBC(field="solute", value = 0, surfaces = [10])
    ]

    model.T = F.Temperature(973)

    derived_quantities = F.DerivedQuantities(
        [F.SurfaceFluxCylindrical(field = "solute", surface=9), 
         F.SurfaceFluxCylindrical(field = "solute", surface=10),
         F.SurfaceFluxCylindrical(field = "solute", surface=8)],
        filename=f"{folder}/derived_quantities.csv",
    )

    model.exports = [F.XDMFExport("solute", folder = folder),
                     derived_quantities
                    ]

    model.settings = F.Settings(
        absolute_tolerance=1e10,
        relative_tolerance=1e-10,
        chemical_pot=True,
        maximum_iterations = 100,
        final_time = 30 * 3600
    )

    model.dt = F.Stepsize(initial_value = 10, stepsize_change_ratio = 1.1)


    model.log_level = 20

    return model

if __name__ == "__main__":


    folder = '3D_cad'
    model = make_model(folder)

    model.initialise()
    model.run()

    data = np.genfromtxt("3D_cad/derived_quantities.csv", delimiter = ',', names=True)
    ts = data['ts']
    flux = data['solute_flux_surface_9']

    plt.plot(ts, flux)
    plt.show()

The transient flux coming from the side wall gives a somewhat-maybe-believable result:

But the flux from the top of the salt gives nonsense:

And even when making a much finer mesh the results become more suspect:
Side flux:

and top flux:

but solving the steady-state solution seems reasonable:

Again would appreciate any input, I’ve been experimenting with lots of different mesh refinements and playing around with the tolerances but nothing comes close to what I would expect the actual results to look like. Thanks.

[remidm]

This is the XDMF solution in transient that I get.

The problem is that your mesh is not in metres (probably in millimetres). See below:

In SALOME, you need to set the dimensions of your geometry to metres otherwise, FESTIM will understand that there is 2 metres of nickel to go through

What happened here is that with the stepsize that you gave, the mesh wasn’t refined enough. But if you increase the stepsize dramatically and also increase the final time, you will see a correct diffusion behaviour.

Consider scaling your geometry in SALOME.