Symmetry boundary condition

User Support
1

Hi all,

I am working with a model that includes three different materials: tungsten, Eurofer, and FLiBe. The boundary conditions are shown in the picture below:

image
image1918×913 393 KB

I am trying to simulate an advection–diffusion problem in which the FLiBe moves with a parabolic velocity profile. However, when I plot the mobile concentration, I notice that, instead of exiting through the top of the FLiBe region as expected, the concentration appears to leave through the Eurofer–FLiBe interface.
image
image999×1125 96.5 KB

I suspect that this issue arises because no boundary condition has been imposed at the top of the FLiBe region. By default, the software applies a Neumann boundary condition with zero flux, which likely prevents the concentration from exiting there.

To address this, I tried imposing a Neumann boundary condition using the known concentration in the FLiBe (previous plot 2.832e+19 [H/m^3]) and the average velocity (1.33 [m/s]). However, the resulting solution is affected by numerical errors:

Neum_con=F.FluxBC(surfaces=12, value=2.832e+19*1.33, field=0) #[H/m^2/s]

image
image998×1125 100 KB

I am not sure how to correctly impose a symmetry boundary condition so that the concentration can leave from the top boundary. Could you please help me with this?
I’m also attaching the code in case someone would like to reproduce the issue.

Thank you in advance for your help.

Best regards.

mesh_domains.h5 (326.7 KB)
mesh_domains.xdmf (585 Bytes)
mesh_boundaries.xdmf (588 Bytes)
mesh_boundaries.h5 (195.1 KB)

import numpy as np
import meshio
import matplotlib.pyplot as plt
import festim as F
import os
import copy
from fenics import plot

my_model = F.Simulation()

my_model.mesh = F.MeshFromXDMF(volume_file=“…/mesh_domains.xdmf”, boundary_file=“…/mesh_boundaries.xdmf”)

tungsten = F.Material(
id=8,
D_0=4.1e-7/(np.sqrt(3)),
E_D=0.39,
S_0=1.35e22,
E_S=0.29,

)

euro_fer = F.Material(
id=9,
D_0=4.57e-7/(np.sqrt(3)),
E_D=0.23,
S_0=1.35e22,
E_S=0.16,
)

FLiBe = F.Material(
id=10,
D_0=9.3e-7,
E_D=42.5/96.48,
S_0=(7.9e-2)*(6.022e23),
E_S=35/96.48,
)

my_model.materials = [tungsten, euro_fer, FLiBe]
from festim import x
import sympy as sp

m1 = (1308.1 - 1337.0) / 0.01 # tra 0 e 0.01
m2 = (1273.0 - 1308.1) / 0.03 # tra 0.01 e 0.04

T_expr = sp.Piecewise(
(1337.0 + m1 * x, x <= 0.01),
(1308.1 + m2 * (x - 0.01), (x > 0.01) & (x <= 0.04)),
(1273.0, True)
)

my_model.T = F.Temperature(T_expr)

implatantion_flux=F.ImplantationDirichlet(surfaces=6,phi=7.15e12, R_p=9.52e-10, D_0=tungsten.D_0, E_D=tungsten.E_D)

Neum_con=F.FluxBC(surfaces=12, value=2.832e+19*1.33, field=0)

Diriclet_con=F.DirichletBC(surfaces=7, value=4.6e-5*(6.022e23), field=0)

my_model.boundary_conditions=[implatantion_flux, Diriclet_con, Neum_con]

my_model.sources=[F.Source(volume=10, field=0, value=1.1e19)]

my_model.settings=F.Settings(
absolute_tolerance=1e-6,
relative_tolerance=1e-6,
transient=False,
#final_time=1000
)
“”"
my_model.dt=F.Stepsize(
initial_value=1e-4,
stepsize_change_ratio=1.2,
max_stepsize=1,
dt_min=1e-9,
milestones=[10, 100, 500, 1000]
)
“”"
my_model.initialise()
my_model.run()
hydrogen_concentration = my_model.h_transport_problem.mobile.solution
cs=plot(hydrogen_concentration, title=" H transport mobile pure diffusion")
plt.colorbar(cs)
plt.show()

from fenics import interpolate, UserExpression, VectorFunctionSpace

functionspace = VectorFunctionSpace(my_model.mesh.mesh, “CG”, 1)

class ParabolicVelocity(UserExpression):
def eval_cell(self, value, x, ufc_cell):
# x[0] = coordinata x, x[1] = y
if x[0] < 0.04 or x[0] > 0.06:
value[0] = 0.0 # componente x della velocità
value[1] = 0.0 # componente y (verticale)
else:
vmax = 2.0
xc = 0.05
r = 0.01 # raggio del profilo parabolico
vy = vmax * (1 - ((x[0] - xc) / r) ** 2)
value[0] = 0.0 # flusso verticale, nessuna velocità orizzontale
value[1] = vy # componente verticale
def value_shape(self):
return (2,) # perché è un vettore (vx, vy)

velocity = ParabolicVelocity(degree=2)

velocity = interpolate(velocity, functionspace)
CS = plot(velocity)
plt.colorbar(CS)
plt.show()
from dolfin import XDMFFile

folder = “…”

filepath = os.path.join(folder, “velocity_field.xdmf”)

xdmf = XDMFFile(filepath)
xdmf.write(velocity)
xdmf.close()

from fenics import inner, dot, grad
my_derived_quantities = F.DerivedQuantities(
[ F.TotalVolume(field=“solute”, volume=10),
F.AverageVolume(field=“solute”, volume=10),
F.PointValue(field=“solute”, x=[0.05, 0.1]),])

my_model.exports = [ my_derived_quantities]

my_model.initialise() # reinitialisation is needed

hydrogen_concentration = my_model.h_transport_problem.mobile.solution
test_function_solute = my_model.h_transport_problem.mobile.test_function

advection_term = inner(dot(grad(hydrogen_concentration), velocity), test_function_solute) * my_model.mesh.dx(10)

my_model.h_transport_problem.F += advection_term
print(my_model.h_transport_problem.F)

my_model.run()
cp=plot(hydrogen_concentration,title=" H transport mobile advection diffusion")
plt.colorbar(cp)
#plot(velocity)
plt.show()

2

Cool model!

In this problem, the right BC at the top/outlet of the FLiBe volume is a “do-nothing” BC so putting the default one is correct. An homogeneous Neumann BC enforces the gradient of concentration to be zero, that doesn’t mean that H isotopes cannot leave, the total flux is the diffusion flux + advective flux.

The first screenshot you show makes sense to me. To convince yourself, you can try and calculate the inlet flux vs outlet advective fluxes with:

J = \int_{\Gamma} c \ v \cdot \mathbf{n} \ ds

You can try and calculate with a back of the envelope calculation or by writing your own custom DerivedQuantity class similar to what is done in this tutorial

3

Also, for this sort of application, I strongly recommend using FESTIM2 which comes with nice features for multi-domain and multi physics problems.