Modeling contamination diffusion in a draining container, part 2

Question

Part 1, but I'll repeat here. This time we'll move the top boundary.

I have a container of water with initial volume $V_0$ and height $L$, open at the top which receives a constant mass flux $\phi$. A sampling line is placed in the container all the way to the bottom which drains the container at some rate $f$. This means the water level drops at a speed $v = fL/V_0$. Thus, in terms of dimensionless time $T = Dt/L^2$ and position $X=x/L$,

$$ \frac{\partial C}{\partial T}=\frac{\partial^2 C}{\partial X^2}+P\frac{\partial C}{\partial X},\quad P=\frac{vL}{D} $$

where I use $+$ for the advection term because $X=0$ will be the bottom of the container, $X=L$ will be the top, and the velocity is in the $-X$ direction. Because the top boundary is moving as the water drains, we can change coordinates to keep the spatial grid fixed instead. I'm adapting the procedure from this paper. We have

$$ X \in [0,\delta(T)],\quad\delta(T) = 1-PT\\[1em] Y = \frac{X}{\delta},\quad Y\in[0,1] $$

Now we transform the original PDE for $C(X,T)\to\tilde C(Y,T)$ with

$$ \frac{\partial C}{\partial X}=\frac{1}{\delta}\frac{\partial \tilde C}{\partial Y}\\[1em] \frac{\partial^2 C}{\partial X^2}=\frac{1}{\delta^2}\frac{\partial^2 \tilde C}{\partial Y^2}\\[1em] \frac{\partial C}{\partial T}=\frac{\partial \tilde C}{\partial T} + \frac{\partial \tilde C}{\partial Y}\frac{\partial Y}{\partial T} = \frac{\partial \tilde C}{\partial T} + \frac{\partial \tilde C}{\partial Y}(-\frac{Y}{\delta})(-P) $$

which we substitute in to get,

$$ \frac{\partial \tilde C}{\partial T} = \frac{1}{\delta^2}\frac{\partial^2 \tilde C}{\partial Y^2} + \frac{P(1-Y)}{\delta}\frac{\partial \tilde C}{\partial Y} $$

Now my question is how to set the boundary conditions. In Part 1, user Rigel confirmed that the bottom boundary should be

$$\frac{\partial C}{\partial X}\Bigg\rvert_{X=0} = 0$$

which lets concentration flow out but not diffuse out. The top boundary I still don't understand. There is a constant flux coming in, which was handled in Part 1 by using the flow like $vC(L,T) = \phi$, but there is no longer any velocity at the top boundary with these coordinates as $P(1-Y)\to 0$ as $Y\to 1$. If there's no velocity, and there can't be diffusion across the top, do we have to use a source term to introduce the mass flux?

Answer 1

The condition you need for the top boundary is a Neumann condition, because the contaminant is entering into the domain at a prescribed rate. The condition is obtained be equating the diffusive flow of contaminant at the boundary with the fixed rate $\phi$:

$$D \frac{\partial C}{\partial x} = \phi$$

Normalizing and transforming it into your reference frame, it becomes (I hope I got the normalization right!)

$$ \frac{\partial \tilde C}{\partial Y} \frac{1}{\delta} = \frac{L}{D} \phi $$

Strictly speaking the water surface is escaping the contaminant. If the contaminant is inflowing in the bottle with velocity $v_c$ and linear density $\rho=\phi/v_c$, such that $v_c \rho = \phi$, then we have that the actual flow of contaminant entering the water is

$$\phi' = (v_c - v) \rho$$

However, if $v_c \gg v$ then $\phi' \approx \phi$ and your top boundary condition writes as above.

Answer 2

In case this helps anyone in the future, I wanted to finish off the problem by showing how I actually implemented it numerically. Forget Newton's method; I just solved the matrix equation for backwards Euler directly, which works great for my parameters anyway ($P \sim0.5$ to $1.5$).

The grid is discretized in $N$ points from $i = 0$ to $N-1$, $\Delta Y = 1/(N-1)$. The PDE on the interior points becomes (drop the tilde on $C$),

$$ \frac{C_{i}^{j+1}-C_i^j}{\Delta T} = \frac{1}{\delta^2\Delta Y^2}\left(C_{i-1}^{j+1}-2C_i^{j+1}+C_{i+1}^{j+1}\right) + \frac{P(1-i\Delta Y)}{\delta\Delta Y}\left(C_{i+1}^{j+1} - C_i^{j+1}\right) $$

and using centered differences at the boundaries, the Neumann BCs give us at $Y=0$:

$$\frac{C_{0}^{j+1}-C_0^j}{\Delta T} = \frac{1}{\delta^2\Delta Y^2}\left(2C_{1}^{j+1}-2C_0^{j+1}\right) + \frac{P}{\delta\Delta Y}\left(C_{1}^{j+1} - C_0^{j+1}\right) $$

and at $Y=1$ (the convection term dies):

$$\frac{C_{N-1}^{j+1}-C_{N-1}^j}{\Delta T} = \frac{1}{\delta^2\Delta Y^2}\left(2C_{N-2}^{j+1}-2C_{N-1}^{j+1}+\frac{2\Delta Y \delta L}{D}\phi\right) $$

The implicit equation tells us then that

$$ \mathbf{AC}^{j+1} = \mathbf{C}^j+\mathbf{b} $$

$\mathbf{A}$ is a tridiagonal matrix where for the interior points:

\begin{align} &\mathbf{A}_{i,i} = 1+2k_1 + (1-i\Delta Y)k_2\\ &\mathbf{A}_{i,i-1} = -k_1\\ &\mathbf{A}_{i,i+1} = -\left(k_1 + (1-i\Delta Y)k_2\right)\\ \end{align}

and for the first and last row:

\begin{align} &\mathbf{A}_{0,0} = 1+2k_1 + k_2\\ &\mathbf{A}_{0,1} = -(2k_1 + k_2)\\ &\mathbf{A}_{N-1,N-1} = 1+2k_1\\ &\mathbf{A}_{N-1,N-2} = -2k_1 \end{align}

where

$$ k_1 = \frac{\Delta T}{\delta^2 \Delta Y^2},\quad k_2 = \frac{P\Delta T}{\delta\Delta Y} $$

And $\mathbf{b}$ is all zeros except for the righthand BC:

$$ \mathbf{b}_{N-1} = 2k_1\delta\Delta Y\frac{L\phi}{D} $$

For my parameters, $\Delta T = $ 5E-4 works fine and I might could go higher. There may be a slight decay in the incoming mass, so double-check what you're running.