I am looking for the analytical solution of 2D Advection-diffusion equation on a rectangular domain $\Omega= (0, a)\times(0, b)$ considering: \begin{equation} \begin{split} & C_t +\lambda C + \nabla\cdot\mathbf{u}C - \nabla\cdot(\mathbf{D} \nabla C) = 0, \text{ for } (x, y, t) \in \Omega \times \mathbb{R^*_+} , \\ & C(x,y,t=0) = \rho(x,y) = f(x)g(y), \text{ for } (x,y) \in \Omega,\\ & f(x) =\left\{ \begin{array}{cl} 1, &x \in [x_0,x_1]\\ 0, & \text{elsewhere}\\ \end{array} \right. g(y) = \left\{ \begin{array}{cl} 1, &y \in [y_0,y_1]\\ 0, & \text{elsewhere}\\ \end{array} \right. \\ & C(x=0,y,t) = C(x=a,y,t) = 0, t \in \mathbb{R^*_+},\\ & \partial_y C(x,y=0,t) = \partial_y C(x,y=b,t) = 0, t \in \mathbb{R^*_+},\\ & \mathbf{u}=(u_x=0,u_y=cst>0) \label{problem} \end{split} \end{equation}
The objective of my post is to verify that the solution calculated with a numerical solver (OpenFOAM...) comply with the analytical solution.
I have successfully managed to find the solution for pure Dirichlet BC and verified the compliance between the numerical solver, but I have some "doubts" when I introduce the Neumann BC $\partial_y C(x,y=0,t) = \partial_y C(x,y=b,t) = 0$.
My very first question is: is the problem well-posed? I had a look on related topic and I found this. However my problem is a bit different than the one described and I don't really know if the answer is applicable here. My concern is that since $C(x,y=0,t) = C(x,y=b,t) \neq 0$, intuitively, we expect "mass injection" at the $C(x,y=0,t)$ for $t \in \mathbb{R^*_+}$. Having tried to solve the problem with OpenFOAM, I get quickly a "diverging solution" where the mass is increasing seemingly indefinitely. Having tried with another solver the solution does not diverge, despite the "mass injection". I would like to decide between the solvers.
My second question holds only if the problem is well-posed. In order to find the solution I used the classical separation of variables method $C(x,y,t)=X(x)Y(y)T(t)$. When I want to solve $Y(y)$, I consider that the form of the solution is: \begin{equation} Y(y)=\exp{(\gamma y)}\big(A\cos{(\eta y)}+B\sin{(\eta y)}\big), \end{equation} where $\gamma\pm i\eta$ are the complex solutions of the corresponding characteristic polynomial: \begin{equation} \begin{split} \gamma\pm i\eta & =\frac{-u_x}{-2D_x}\pm i\frac{\sqrt{|\Delta|}}{-2D_x}, \\ & =\frac{u_x}{2D_x}\pm i\frac{\sqrt{|\Delta|}}{2D_x}. \end{split} \end{equation} With the the previously defined BC, I have to solve the following system: \begin{equation} \begin{split} \partial_yY(0)=\exp{(\gamma 0)}\big((\gamma A+\eta B)\cos{(\eta 0)}+(\gamma B-\eta A)\sin{(\eta 0)}\big) & = 0,\\ \partial_yY(b)=\exp{(\gamma b)}\big((\gamma A+\eta B)\cos{(\eta b)}+(\gamma B-\eta A)\sin{(\eta b)}\big) & = 0. \end{split} \end{equation} The solutions of the previous system are: \begin{equation} \begin{split} A = -\frac{\eta B}{\gamma },\\ B\bigg(\gamma+\frac{\eta^2}{\gamma}\bigg)\sin{(\eta b)} & = 0, \gamma \neq 0. \end{split} \end{equation} In order to satisfy the second equation of the previous system I choose (considering a Fourier series to come): \begin{equation} \begin{split} \eta &= \eta(m) = \eta_m = \frac{m\pi}{b},\\ B & = B(m)=B_m,~m\in\mathbb{N^*}. \end{split} \end{equation} The form of the solutions is: \begin{equation} Y_0(y)= B_m\exp{(\gamma y)}\bigg(-\frac{\eta_m}{\gamma }\cos{(\eta_m y)}+\sin{(\eta_m y)}\bigg),~m\in\mathbb{N^*}, \gamma \neq 0. \end{equation}
I have a question at this stage, unlike for the Dirichlet BC solution, I notice that $\gamma\neq0$ which means that the solution does not stand anymore if $u_y=0$. How should I interpret this considering the concern express in my first question?
Pursuing with the normalization I have to calculate the scalar product of two eigenfunctions which leads formally to: \begin{equation} \begin{split} N^2_{jk} & = \int_0^b \bigg(-\frac{j\pi}{b\gamma}\cos{(\frac{j\pi y}{b})}+\sin{(\frac{j\pi y}{b})}\bigg) \bigg(-\frac{k\pi}{b\gamma }\cos{(\frac{k\pi y}{b})}+\sin{(\frac{k\pi y}{b})}\bigg) dy,\\ N^2_{jk} & = \bigg(\frac{jk\pi^2}{2b\gamma^2} + \frac{b}{2}\bigg)\delta_{jk},\\ N_{j=k} & = \sqrt{\frac{j^2\pi^2}{2b\gamma^2} + \frac{b}{2}}. \end{split} \end{equation} I notice that norm won't be constant with regard to $m$ (or $j,k$) for the Fourier series like with the Dirichlet BC. Again, How should I interpret this?
$B_m$ Fourier series coefficients are determined through: \begin{equation} \begin{split} B_m&=\frac{1}{N_{j=k=m}}\int_{0}^b \exp{(-\gamma y)}g(y)\bigg(-\frac{\eta_m}{\gamma }\cos{(\eta_m y)}+\sin{(\eta_m y)}\bigg)d y,\\ B_m & = \bigg(\frac{m^2\pi^2}{2b\gamma^2} + \frac{b}{2}\bigg)^{-\frac{1}{2}}\bigg[ -\frac{1}{\gamma}\exp(-\gamma y)\sin{(\eta_m y)}\bigg]_{y_0}^{y_1}, \end{split} \end{equation}
Finally, when I try to calculate the solution it doesn't seem to give consistent results.
Any help would be greatly appreciated.
Best regards,
In case this may help some people, after trials and errors I finally figured out that the normalization factor $N^2_m$ had to be applied...
If I trust the calculation, I guess I can also conclude that the problem is well-posed...
Best regards,