Understanding units for TotalVolume

User Support
1

Hi all, I’m investigating T inventory in a pipe. I originally used a 3D pipe and all was working fine until I played around with surface parameters and it began struggling to converge. I tried adjusting tolerances but it was very finickity so I switched to a 2D pipe slice, with a finer mesh. Obviously my results are now different as it’s finding the inventory in a 2D area.

Before, for the 3D pipe (25mm ID, 35mm OD, L=100mm), I assumed the result was the total number of mols in the volume, so was multiplying the output by 10 to get mols per meter of pipe.

For the 2D simulation, how do I interpret the values I obtain? And how do I normalise to mol per meter of pipe? I took a look here but am still a little confused.

Many thanks!

Stefan

Edit: for clarity, here are my results:

image
image1165×650 42.6 KB

image
image1791×855 75.1 KB

You can see that the 3D simulation is giving me a smaller output for TotalVolume. I assumed this was in the unites mol, so I multiplied by 10 to get mol per meter of pipe (as my mesh was a 3D pipe of 0.1m).

As for the 2D results, I’m not sure what the units are, or how to convert them to mol per meter of pipe.

2

Hi! So just to clarify: you are running a 2D simulation of a pipe and computing TotalVolume on it.

Is it the radial cross section or an axial cross section? Since you’re asking about “per meter of pipe” I’m assuming it’s an axial cross section like this

image
image888×719 26.7 KB

If you calculate TotalVolume on this then the units should already be per meter of pipe.

As you rightfully pointed out, TotalVolume integrates \int c dA where c is still in mol/m3 and dA is a surface element in m2. The quantity is therefore in mol/m.

Let me know if there is something I’m not getting.

Remi

3

Thanks for your response - yes, the 2D simulation is an axial slice. So it seems TotalVolume for 2D is mol per meter of pipe, and TotalVolume for 3D is mol in the volume provided (i.e. 0.1m length)?

If this is the case, multiplying my results from the 3D simulation should give the same results as the 2D simulation, however this is currently not the case. I wonder why this is - it could be the difference in mesh? It’s quite a large difference between the two!

4

This is correct.

It is a bit hard to assist you without seeing your code. Could you post a MWE (minimal working example) reproducing what you are seeing?

5

@srichr here’s a small example showing that things should work:

from mpi4py import MPI
import gmsh, numpy as np, festim as F
from dolfinx.io.gmsh import model_to_mesh

R_in, R_out, L = 12.5e-3, 17.5e-3, 100e-3  # m
source, final_time = 1e18, 100            # H/m3/s, s

def pipe_mesh(dim, h=1e-3):
    """2D: annulus cross-section. 3D: the same annulus extruded over L."""
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0)
    outer = gmsh.model.occ.addDisk(0, 0, 0, R_out, R_out)
    inner = gmsh.model.occ.addDisk(0, 0, 0, R_in, R_in)
    entities = gmsh.model.occ.cut([(2, outer)], [(2, inner)])[0]
    if dim == 3:
        entities = [e for e in gmsh.model.occ.extrude(entities, 0, 0, L) if e[0] == 3]
    gmsh.model.occ.synchronize()
    gmsh.model.addPhysicalGroup(dim, [tag for _, tag in entities], 1)
    gmsh.option.setNumber("Mesh.MeshSizeMax", h)
    gmsh.model.mesh.generate(dim)
    mesh = model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=dim).mesh
    gmsh.finalize()
    return mesh

def run(dim):
    my_model = F.HydrogenTransportProblem()
    my_model.mesh = F.Mesh(pipe_mesh(dim))
    vol = F.VolumeSubdomain(id=1, material=F.Material(D_0=1e-7, E_D=0.2))
    H = F.Species("H")
    my_model.subdomains, my_model.species = [vol], [H]
    my_model.temperature = 600
    my_model.sources = [F.ParticleSource(value=source, volume=vol, species=H)]
    total = F.TotalVolume(field=H, volume=vol)
    my_model.exports = [total]
    my_model.settings = F.Settings(atol=1e8, rtol=1e-10, final_time=final_time)
    my_model.settings.stepsize = F.Stepsize(final_time / 2)
    my_model.initialise()
    my_model.run()
    return total.value

total_2d, total_3d = run(2), run(3)
inventory = source * final_time * np.pi * (R_out**2 - R_in**2)  # H per m of pipe
print(f"2D: {total_2d:.4e} H/m   (expected {inventory:.4e})")
print(f"3D: {total_3d:.4e} H     (expected {inventory * L:.4e})")
print(f"3D / 2D = {total_3d / total_2d:.4f} m   (pipe length = {L} m)")

Produces

2D: 4.7123e+16 H/m   (expected 4.7124e+16)
3D: 4.7124e+15 H     (expected 4.7124e+15)
3D / 2D = 0.1000 m   (pipe length = 0.1 m)
6

@srichr have you been able to fix your issue?

7

Hi @remidm, thank you for your responses. Apologies for my delayed response; I’ve been off sick. My colleague thinks it could be due to convergence as a result of me having been using an older version of FESTIM2. I’ve updated to the latest release; let me dig into what’s going on a little then I’ll post a MWE.

Thanks!
Stefan

8

@remidm the mesh was to blame! :man_facepalming: I re-made the 2D mesh and results are matching up once again. Thanks for your help.

Stefan