Hello,
I am new to fenics and FEM simulation in general, and am trying to run a “simple” 2D problem consisting of a solution with some species and 2 electrodes in it. I want to get the concentration profile in the bulk solution as the species gets consumed at an electrode and remaining species start diffusing toward the electrode. Since the species cannot diffuse through the electrodes themselves, I used a dirichlet condition on concentration to set it to 0 within the electrodes area (which are rectangles). At each timestep I enforce this c = 0 condition. I also start with the initial condition that concentration is 0 in the electrodes area and 12 in the bulk solution.
Here’s how I solve my system:
a = u*v*dx + D*dt*inner(grad(u),grad(v))*dx
L = c_0 * v * dx
u = Function(V)
bc_electrode_left = DirichletBC(V, Constant(0.0), electrode_left)
bc_electrode_right = DirichletBC(V, Constant(0.0), electrode_right)
bcs = [bc_electrode_left,bc_electrode_right]
solve(a == L, u, bcs)`
My problem is that at t = 0 the solution looks like this:
(which is fine to me)
However the next timestep it looks like this:
So some diffusion is happening while I don’t consume anything at the electrodes. I read it could be because I didn’t set a neumann condition on the flux normal to the electrode to make it an impermeable barrier, but even when reading examples/tutorials using such conditions I could not figure out how to do it.
I believe I could do it by excluding the electrodes from the mesh (I read that the borders of the domain are by default impermeable?), but problem is I also need the electrodes to grow during the simulation as matter deposits on it, and I believe it is not possible to change the domain at each timestep?
Therefore I’d like to request help on how to make my electrodes completely impermeable to diffusion, I’d appreciate any pointers you could have
Edit: here’s the minimal code (I tried to get rid of most of the unrelated stuff)
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
import mshr as mshr
import shapely as shp
import math
def point_inside_polygon(point, polygon):
polygon = shp.Polygon(polygon)
return polygon.contains(shp.Point(point[0],point[1]))
# Define the geometry parameters
length_domain = 440.0e-6 # µm
width_domain = 225.0e-6 # µm
length_electrode = 100.0e-6 # µm
width_electrode = 25.0e-6 # µm
gap_width = 240.0e-6 # µm
start_y = 100.0e-6
class ElectrodeBoundaryLeft(SubDomain):
def __init__(self, initial_points):
super().__init__()
self.perimeter_points = initial_points
def inside(self, x, on_boundary):
return point_inside_polygon(x, self.perimeter_points)
class ElectrolyteBoundary(SubDomain):
def __init__(self, left_elec):
super().__init__()
self.left_elec = left_elec
def inside(self, x, on_boundary):
return not (self.left_elec.inside(x, on_boundary))
perimeter_points_left = [(0, start_y), (length_electrode, start_y), (length_electrode, start_y + width_electrode), (0, start_y + width_electrode),(0, start_y)]
electrode_left = ElectrodeBoundaryLeft(perimeter_points_left)
D = 1e-6
electrolyte_boundary = ElectrolyteBoundary(electrode_left)
domain = mshr.Rectangle(Point(0, 0), Point(length_domain, width_domain))
left_elec_polygon = mshr.Polygon([Point(x, y) for x, y in electrode_left.perimeter_points])
domain.set_subdomain(1, left_elec_polygon)
mesh = mshr.generate_mesh(domain, 100)
V = FunctionSpace(mesh, 'P', 1)
# Mark the boundaries
boundaries = MeshFunction('size_t', mesh,1)
boundaries.set_all(0)
electrode_left.mark(boundaries, 1)
initial_c = 12.0e-3
# Define initial condition
class InitialCondition(UserExpression):
def __init__(self, electrolyte_bound):
super().__init__()
self.electrolyte_bound = electrolyte_bound
def eval(self, value, x):
if self.electrolyte_bound.inside(x, True):
value[0] = initial_c
else:
value[0] = 0.0
c_0 = interpolate(InitialCondition(electrolyte_boundary), V)
dt = 1e-3
t = dt
T = 2
while t < T:
u = TrialFunction(V)
v = TestFunction(V)
a = u * v * dx + D*dt*inner(grad(u), grad(v))*dx
L = c_0 * v * dx
u = Function(V)
bc_electrode_left = DirichletBC(V, Constant(0.0), electrode_left)
bcs = [bc_electrode_left]
solve(a == L, u, bcs)
fig = plt.figure()
ploooot = plot(u, title="Concentration")
plt.colorbar(ploooot)
plt.savefig("field/{:04f}_concentration.png".format(t))
vtkfile = File("field/{:04f}_concentration.pvd".format(t))
vtkfile << u
t = t + dt
c_0.assign(u)
plt.close('all')