Hi all,
I’m interested on modelling an electrolytic pickling process. The mathematical description of this problem is the following. Consider the geometry:
I would like to solve:
\nabla \cdot (-k_i\nabla u) = 0 \hspace{1.4cm} \text{in } \Omega_i; \quad i = 1, 2
\hspace{.65cm} - k\nabla u \cdot \nu = g(u) \hspace{.75cm} \text{in } \partial \Omega
\hspace{.75cm} [k\nabla u \cdot \nu] = 0 \hspace{1.38cm} \text{in } \Gamma_I
\hspace{.65cm}- k\nabla u \cdot \nu = h([u]) \quad \text{in } \Gamma_I
With \nu pointing towards \Omega_2 and where both g(u) and h([u]) are nonlinear functions on their respective arguments. In particular I have in mind functions of the form:
g(u) = \begin{cases} i_{a,0}\exp(\alpha_a(U_a - u - E_{a,eq})) && \text{in } \Gamma_a \\ -i_{c,0}\exp(-\alpha_c(U_c - u - E_{c,eq})) && \text{in } \Gamma_c \\ 0 && \text{otherwise} \end{cases}
and
h([u]) = \begin{cases} i_{a,0}^*\exp(\alpha_a^*([u] - E_{a,eq})) && \text{if } [u] \ge E_{a,eq}\\ -i_{c,0}^*\exp(-\alpha_c^*([u] - E_{c,eq})) && \text{if } [u] \le E_{c,eq} \\ 0 && \text{otherwise} \end{cases}
I attempted to solve this problem by employing the (I think) equivalent weak formulation: Find (u_1, u_2, i_a, i_c, i) \in W:
\int_{\Omega_1} k_1\nabla u_1 \cdot \nabla v_1dx + \int_{\Gamma_I} iv_1 ds = 0
\int_{\Omega_2} k_2\nabla u_2 \cdot \nabla v_2dx - \int_{\Gamma_I} iv_2 ds - \int_{\Gamma_a} i_av_2ds - \int_{\Gamma_c} i_cv_2= 0
\int_{\Gamma_I}i\eta ds = \int_{\Gamma_I}h(u_1 - u_2)\eta ds
\int_{\Gamma_a}i_a\eta_ads = \int_{\Gamma_a}g(u_2)\eta_ads
\int_{\Gamma_c}i_c\eta_cds = \int_{\Gamma_c}g(u_2)\eta_cds
for any (v_1, v_2, \eta_a, \eta_c, \eta) \in W
I attempted to implement it by using the Mixed-dimensional branch by ceciledc, with the following code:
from fenics import *
from mshr import *
import logging
logging.getLogger('FFC').setLevel(logging.WARNING)
#Parameters
k_0 = 1e7
k_1 = 10
i_anode_0 = 5e-5
i_cath_0 = 7.9e-4
i_an_b_0 = 5e-10
i_ca_b_0 = 3.2e-3
alpha_a = 13.54
alpha_c = 16.33
alpha_a_b = 28.27
alpha_c_b = 13.14
U = 8.65
E_a_eq = 0.82
E_c_eq = -0.41
#Creation of mesh
background = Rectangle(Point(0,0), Point(0.3,0.07))
steel = Rectangle(Point(0,0), Point(0.3,0.01))
anode = Rectangle(Point(0.025,0.025), Point(0.125,0.035))
cathode = Rectangle(Point(0.175,0.025), Point(0.275,0.035))
electrolyte = background - steel - anode - cathode
domain = steel + electrolyte
domain.set_subdomain(1, electrolyte)
mesh = generate_mesh(domain, 16)
#Define subdomains
class Gamma_a(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and between(x[0], [0.025, 0.125]) and between(x[1], [0.025, 0.035])
class Gamma_c(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and between(x[0], [0.175, 0.275]) and between(x[1], [0.025, 0.035])
class Gamma_I(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.01)
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x[1] <= 0.01
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[1] >= 0.01
#Initializate subdomains instances
gamma_a = Gamma_a()
gamma_c = Gamma_c()
gamma_I = Gamma_I()
omega_0 = Omega_0()
omega_1 = Omega_1()
#Creation of facet and cell markers and marking of mesh elements
facet_marker = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
gamma_a.mark(facet_marker, 1)
gamma_c.mark(facet_marker, 2)
gamma_I.mark(facet_marker, 3)
cell_marker = MeshFunction("size_t", mesh, mesh.topology().dim())
omega_0.mark(cell_marker, 0)
omega_1.mark(cell_marker, 1)
xdmffile = XDMFFile('mesh_function.xdmf')
xdmffile.write(facet_marker)
#Definition of new meshes
mesh_0 = MeshView.create(cell_marker, 0)
mesh_1 = MeshView.create(cell_marker, 1)
mesh_a = MeshView.create(facet_marker, 1)
mesh_c = MeshView.create(facet_marker, 2)
mesh_I = MeshView.create(facet_marker, 3)
#Definition of measures
dx0 = Measure("dx", domain=mesh_0)
dx1 = Measure("dx", domain=mesh_1)
dxa = Measure("dx", domain=mesh_a)
dxc = Measure("dx", domain=mesh_c)
dxI = Measure("dx", domain=mesh_I)
#Definition of function space, trial and test functions
V0 = FunctionSpace(mesh_0, "CG", 2)
V1 = FunctionSpace(mesh_1, "CG", 2)
V2 = FunctionSpace(mesh_a, "CG", 2)
V3 = FunctionSpace(mesh_c, "CG", 2)
V4 = FunctionSpace(mesh_I, "CG", 2)
W = MixedFunctionSpace(V0, V1, V2, V3, V4)
u = Function(W)
(u0, u1, ia, ic, i) = u.split()
(v0, v1, ea, ec, e) = TestFunctions(W)
#Definition of tafel law
def i_tafel(i0, alpha, eta):
return i0*exp(alpha*eta)
#Preliminary expressions
eta_O2 = u0 - u1 - E_a_eq
eta_H2 = u0 - u1 - E_c_eq
H_O2 = conditional(le(eta_O2, 0), 0, 1)
H_H2 = conditional(ge(eta_H2, 0), 0, 1)
i_anode = i_tafel(i_anode_0, alpha_a, U*.5 - u1 - E_a_eq)
i_cathode = i_tafel(-i_cath_0, -alpha_c, -U*.5 - u1 - E_c_eq)
i_strip_a = i_tafel(i_an_b_0, alpha_a_b, eta_O2)
i_strip_c = i_tafel(-i_ca_b_0, -alpha_c_b, eta_H2)
i_interface = i_strip_a*H_O2 + i_strip_c*H_H2
#Variational formulation
int_Omega = inner(k_0*grad(u0),grad(v0))*dx0 + inner(k_1*grad(u1),grad(v1))*dx1
int_anode = ia*(ea - v1)*dxa - i_anode*ea*dxa
int_cathode = ic*(ec - v1)*dxc - i_cathode*ec*dxc
int_interface = i*(v0 - v1 + e)*dxI - i_interface*e*dxI
Res = int_Omega + int_anode + int_cathode + int_interface
#Solution of the system
solve(Res == 0, u)
(u0, u1, ia, ic, i) = u.split()
But I get the following error message:
Error: Unable to evaluate function at point. *** Reason: The point is not inside the domain. Consider calling "Function::set_allow_extrapolation(true)" on this Function to allow extrapolation.
Any idea what the problem might be? I don’t see why there might be an evaluation of a function outside its domain.
P.S. I’m using the docker container of the mixed-dimensional branch:latest in Windows 10.
