2d simulation breaks when introducing traps

User Support
1

Hi,

I am trying to model 2d steel square inside another piece of steel, with traps in the inner piece of steel only. Initially, hydrogen is in inner steel. The code below works ok, but as soon as I uncomment the trap lines, something strange happens: at time = 0, nonzero concentration of solute is at the interface between steels only!

Could you please explain what is wrong with my code?

Best wishes,

Mikhail

from fenics import *

import matplotlib.pyplot as plt

import festim as F

from festim import x, y, z

import numpy as np

import sympy

import math

# GEOMETRY AND MESH

width = 1

height = 1

border_thickness = 0.1

nx, ny = 21, 21

mesh = RectangleMesh(Point(0, 0), Point(width, height), nx, ny)

# MARK VOLUMES

volume_markers = MeshFunction(“size_t”, mesh, mesh.topology().dim())

volume_markers.set_all(2)

plot(volume_markers, title=“Volume Markers (1=steel1, 2=steel2)”)

plt.show()

inner_metal = CompiledSubDomain(

“x[0] > x_min && x[0] < x_max && x[1] > y_min && x[1] < y_max”,

x_min=border_thickness,

x_max=width - border_thickness,

y_min=border_thickness,

y_max=height - border_thickness

)

inner_metal.mark(volume_markers, 1)

model = F.Simulation()

model.mesh = F.Mesh(mesh=mesh, volume_markers=volume_markers)

# MATERIALS

steel1 = F.Material(

id = 1,

D_0 = 1e-16, #m2/s

E_D = 0, #eV

)

steel2 = F.Material(

id = 2,

D_0 = 1E-15, #m2/s

E_D = 0, #eV

)

model.materials = [steel1, steel2]

# TEMPERATURE

model.T = 300

# INITIAL CONDITIONS

model.initial_conditions = [

F.InitialCondition(

    field="solute",

    value=1e26 \* sympy.Piecewise(

        (1, (x > border_thickness) & (x < width - border_thickness) &

            (y > border_thickness) & (y < height - border_thickness)),

        (0, True)

    )

)

]

# SURFACE MARKERS

surface_markers = MeshFunction(“size_t”, mesh, mesh.topology().dim() - 1)

surface_markers.set_all(0)

steel1_boundary = CompiledSubDomain(

“near(x[0], x_min) || near(x[0], x_max) || near(x[1], y_min) || near(x[1], y_max)”,

x_min=border_thickness,

x_max=width - border_thickness,

y_min=border_thickness,

y_max=height - border_thickness

)

steel2_boundary = CompiledSubDomain(

“on_boundary”,

)

steel1_boundary.mark(surface_markers, 1)

steel2_boundary.mark(surface_markers, 2)

model.mesh.surface_markers = surface_markers

SS316L_density = 8.55e28

#TRAPS

#SS316L_trap = F.Trap(

# k_0 = 3.82e-7 / ((5.66e-10)**2 * 6 * SS316L_density),

# E_k = 0.16,

# p_0 = 1e13,

# E_p = 0.5,

# density = 6.08e25,

# materials = [steel1]

#)

#model.traps = [SS316L_trap]

# BOUNDARY CONDITIONS

model.boundary_conditions = [

F.DirichletBC(field = 0,

              value = 0,

              surfaces = \[2\]),

]

#EXPORTS

results_folder = “TWO_DIMENSION_PERMEATION”

model.exports = [

F.XDMFExport(

    field="solute",

    filename=results_folder + "/hydrogen_concentration.xdmf",

    checkpoint=False,

),

]

# SETTINGS

model.settings = F.Settings(

absolute_tolerance = 1e16, 

relative_tolerance = 1e-8,

final_time= 1e3, 

)

# TIMESTEP

model.dt = F.Stepsize(

initial_value=1,

stepsize_change_ratio=1.2, 

dt_min=1e-5,

max_stepsize=1e3,

)

model.initialise()

model.run()

2

P.S. This is FESTIM1

3

Hi @miklavrent could you provide some screenshots to illustrate your issue? The code block isn’t formatted properly so it’s not easy to copy paste

4

Hi Remi,

Thank you for your reply! It looks like there are several unexpected backslashes in the code, perhaps it’ll work if you delete them? Anyway, I attach two screenshots showing initial mobile distribution with and without the trap part.

Best wishes,

Mikhail

with_trap
with_trap1969Ă—957 246 KB

without_trap
without_trap1969Ă—957 163 KB

5

If you look at the colour bar, the first screenshot shows a negative value in the inner square. The red strip is an “overshoot”. This seems like a typical Gibb’s phenomenon issue related to the mesh.

Basically you are trying to approximate very sharp gradients with too few cells. Things to consider:

  • refining the mesh
  • use discontinuous finite elements (DG) for the trap concentration

To get an idea of the refinement level you need, you may want to setup a 1D simulation an play with the refinement there. By the way, are you sure your problem isn’t a 1D problem?

6

Thank you Remi, I’ll try to refine the mesh and/or use DG. Indeed this looks like a 1D problem, but it’s because I simplified it to bring as close to MWE as possible.

1 Like