# Modeling contamination diffusion in a draining container, part 2

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?

• Do you know the velocity ($v_c$) or the linear density ($\rho$) that this inflowing contaminant has? I understand that the contaminant flux is $\phi$ (=contaminant entering over time) but this is the product $\rho v_c$. The top bc depends on them... unless the velocity of the contaminant is much higher than that of water surface. Sep 15 at 16:17
• The velocity is probably much higher; it's blowing down into the water surface. These are extremely fine particles (cleanroom environment) that will just crash into the water surface and stay there, not sink, so would the velocity be discontinuous, at least? But could we just make up a number for $v_c$? Sep 15 at 20:46

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.

• Thank you. At first I thought this wasn't right because it didn't seem to conserve mass correctly (before any concentration reaches the drain, the amount of mass added each timestep should equal the flux), but from experimenting it seems I might have to make $\Delta T$ very small, like 1e-6. Sep 18 at 5:36
• Consider that if you are using an explicit method you are subject to two necessary stability conditions, one for the advective operator and one for the diffusive operator. For the advection you have $\Delta T < \delta \Delta Y / [P (1-Y)]$, and for diffusion $\Delta T < \delta^2 \Delta Y^2$. I don't know if 1e-6 is well below these threshoulds, but you should take them into account anyway. Sep 18 at 14:23
• I'm using backwards Euler, with imaginary points on each end of the domain to control the Neumann BCs and Newton's method to solve each timestep. But I see the same thing: smaller $\Delta Y$ forces me to use smaller $\Delta T$. This is perhaps a separate question, so I can mark this closed. Sep 19 at 2:28

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.