Implementing no-flux boundary condition reaction-diffusion PDE

Question

I'm having trouble figuring out how to implement boundary conditions for this problem: \begin{align} \frac{\partial n}{\partial t} &= D_n\nabla^2n - \nabla\cdot\left(\frac{\chi}{1+\alpha c}n\nabla c\right) - \nabla\cdot (pn\nabla f) \\ \frac{\partial f}{\partial t} &= \beta n-\gamma nf \\ \frac{\partial c}{\partial t} &= -\eta nc \end{align} Where the physical domain is a simple square, $\Omega = [0,100]\times[0,100]$ and the time domain is $T = [0,t]$. The paper I'm working from says the no-flux boundary condition is: $$\left[D_n\nabla n - \left(\frac{\chi}{1+\alpha c}n\nabla c\right) - (pn\nabla f)\right]\cdot \hat{N} = 0$$ Where $\hat{N}$ is the outward unit normal vector.

What I've done: I've set up the central difference versions of the first 3 equations (Overleaf PDF here), but I'm not sure how to set up $\frac{\partial n}{\partial t}$ on the boundary using the BC above. I've tried converting the BC into a finite difference version and using it to cancel some terms, but I'm not sure if this is the right way to go about this. Any advice is appreciated!

Answer 1

A finite volume viewpoint would be very covenient for your problem. It is not too different from finite difference. Sticking to just 1-d, your first cell would be the interval $[0,h]$ and you store the solution at the center of this cell, say $n_0$. Then the semi-discrete scheme for this cell would be $$ \frac{d n_0}{dt} + \frac{F_{1/2} - F_{-1/2}}{h} + ... = 0 $$ From your boundary condition $F_{-1/2} = 0$. The other fluxes can be computed by a finite difference, e.g., $$ F_{1/2} = D_n \frac{n_1 - n_0}{h} - \frac{\chi n_{1/2}}{1+\alpha c_{1/2}} \frac{c_1 - c_0}{h} - p n_{1/2} \frac{f_1 - f_0}{h} $$ where $n_{1/2}=(n_0 + n_1)/2$, etc.

Answer 2

The equations for $f$ and $c$ do not have any spatial derivatives, and consequently they do not require (and cannot have, in fact) any boundary conditions. Thus, only $n$ has boundary conditions, and these are in your case of Neumann type. It is convenient to rewrite the equation you have as $$ \left(\frac{\chi}{1+\alpha c}n\nabla c\right)\cdot\hat{N} = \left[D_n\nabla n - (pn\nabla f)\right]\cdot \hat{N} $$ to make clear what condition has to be satisfied.

In a finite element context, the term on the left hand side of this appears naturally as a result of integration by parts, and dealing with the boundary condition is therefore a trivial matter of substituting the term by what you have here from the boundary condition. I am not enough of a finite difference person to say how best to do this in that context.