# Robin coupled steady-state diffusion (1st order \*heterogeneous\* chemical reaction) **URL:** \ **Category:** variational formulation\ **Created:** [July 6, 2020, 7:31pm UTC](https://fenicsproject.discourse.group/t/robin-coupled-steady-state-diffusion-1st-order-heterogeneous-chemical-reaction/3709 "2020-07-06T19:31:17Z")\ **Posts on this page:** 1\ **Page:** 1
**Author:** ![genericFENICSUser](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/genericfenicsuser/32/1455_2.png) [@genericFENICSUser](https://fenicsproject.discourse.group/u/genericFENICSUser)\ **Post date:** [July 6, 2020, 7:31pm UTC](https://fenicsproject.discourse.group/t/robin-coupled-steady-state-diffusion-1st-order-heterogeneous-chemical-reaction/3709/1 "2020-07-06T19:31:17Z")
**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 : ```python 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](https://fenicsproject.org/qa/10927/neumann-bc-depending-on-solution?show=11000#a11000) 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: [Link 1](https://fenicsproject.org/qa/10927/neumann-bc-depending-on-solution/?show=10927#q10927) [Link 2](https://fenicsproject.discourse.group/t/segregation-flux-boundary-condition-between-two-domains/376)