Hi all, I’ve been trying to solve a coupled diffusion problem but I’m really unsure about how to describe/implement the coupling. I’d appreciate any guidance.
EDIT: Defined the problem a bit better
…
I have a domain \Omega split into sub-domains, \Omega_1 and \Omega_2, with different constant diffusion coefficients (D_1, D_2) that are separated by some internal boundary \Gamma. I am interested in the steady-state solution to Fick’s 2nd law across the system in two dimensions, given constant Dirichlet BC everywhere on \partial \Omega:
0 = D\left(\frac{\partial^2 c}{\partial x^2} + \frac{\partial^2 c}{\partial x^2}\right)
which I have in the form :
a = D_1 * inner(grad(u), grad(v)) * dx(1) \
+ D_2 * inner(grad(u), grad(v)) * dx(2)
f = Constant(0.0)
L = f * v * dx
The part I’m unsure about is the coupling on the internal boundary. The concentrations in each subdomain at \Gamma are coupled by a constant K. Therefore I need to the following BC
\frac{\partial c_{\Gamma}}{\partial n} = K \left(c_{\Gamma, {\Omega_2}} - c_{\Gamma, {\Omega_1}}\right)
Which is a Robin boundary condition of the form shown in the FENICS Tutorial 4.4.1. That should also imply that the flux out of \Omega_1 is equal to the flux into \Omega_2
{\frac{\partial c}{\partial n}}_{\Gamma, \Omega_1} = {\frac{\partial c}{\partial n}}_{\Gamma, \Omega_2}
The implementation I have been trying based on Link 1 where the system is first implemented as mixed continuous Galerkin space and then projected onto a discontinuous Galerkin space. In that case however, a Neumann boundary condition couples the sub-domains.
I’ve been trying this for a while without success so any guidance is much appreciated!
Thanks!
Additional Info
FYI. The problem represents heterogeneous chemical kinetics at an inter-phase boundary where a chemical species passes between the phases. The kinetic term is represented by the Butler-Volmer model from electrochemistry.
There are several similar examples that I have seen on here: