2D Stationary advection-diffusion eq. as a BVP

[mhsd91]

PROBLEM FORMULATION:
Considering the region [itex] \Omega [/itex] bounded as a square box within [itex] x \in [0,1], y \in [0,1] [/itex]. We wish to solve the 2D, stationary, advection-diffusion equation,

[itex]
0 = D\nabla^2 \rho(x,y) + \vec{V} \cdot \nabla \rho(x,y)
[/itex]

where [itex] D [/itex] is a scalar constant, and [itex] \vec{V} = V_1\hat{e}_x + V_2 \hat{e}_y [/itex] is a constant advection field vector. The problem has the following boundary conditions,

[itex]
\begin{align}
\nabla \rho(1,y) &= \nabla \rho(x,1) =& 0 \\
\rho(0,y) &= \rho(x,0) =& C = \textrm{known constant}
\end{align}
[/itex]

ATTEMPT AT SOLUTION: Method of Separation of Variables

Assuming [itex] \rho(x,y) = X(x) \cdot Y(y) [/itex], then

[itex]
\begin{align}
\nabla \rho &= Y X_{x} \hat{e}_x + X Y_{y} \hat{e}_y \\
\nabla^2 \rho &= Y X_{xx} + X Y_{yy}
\end{align}
[/itex]

where the subscript denotes partial derivation with respect to that particular variable, with the exception of [itex]\hat{e}[/itex] which represents the unit vector in either direction. Decomposing the advection field vector and inserting these results into the original problem, we get..

[itex]
0 = D (Y X_{xx} + X Y_{yy} ) + Y X_{x} v_{x} + X Y_{y} v_{y}
[/itex]

Multiply this eq. with [itex] 1/(XYD) [/itex] and moving all [itex] X [/itex] -expressions to the left hand side result in the separation of [itex] X [/itex] and [itex] Y [/itex] such that they both have to be equal some unknown constant [itex] \lambda [/itex] (we cannot have a small change in [itex] X [/itex], without the corrosponding change in [itex] Y [/itex] and vice versa). Thus,

\begin{align}
\lambda &= -\frac{X_{xx}}{X} - \frac{X_{x}}{X} \frac{V_{1}}{D} \\
\lambda &= +\frac{Y_{yy}}{Y} + \frac{Y_{y}}{Y} \frac{V_{2}}{D}
\end{align}

Which are two independent, linear, second order ODEs, with general solutions
[itex]
\begin{align}
X(x) &= C_1 \exp[-\alpha^{(+)}_{x} x] \\
&+ C_2 \exp[+\alpha^{(-)}_{x} x] \\
& \\
Y(y) &= C_3 \exp[-\alpha^{(+)}_{y} y] \\
&+ C_4 \exp[+\alpha^{(-)}_{y} y]
\end{align}
[/itex]

with

[itex]
\alpha^{(\pm)}_{m} = \frac{1}{2} \left( \sqrt{\frac{v_{m}^2}{D^2} +4\lambda} \pm \frac{v_{m}}{D}\right), \quad m=1 \vee 2
[/itex]

and [itex] C_i, i=1,2,3,4; [/itex] are constant coefficients. So far, so good! However, here is where my issues begin to pile up as I'm unable to sort out the [itex] C [/itex] -values with the boundary conditions.

Along the north and east boundaries, I'm able to write [itex] C_1 = (\textrm{some const expression}) \cdot C_2 [/itex] and similar for [itex] C_3, C_4 [/itex]. However, for the west and south bounds I end up with

[itex]
(C_1 + C_2)Y = C \\

(C_3 + C_4)X = C
[/itex]

Which will only be valid for [itex]C=0, C_1 = -C_2, C_3=-C_4[/itex]. This, evidently, results in the trivial solution [itex]\rho(x,y)=0[/itex], which obviously is not what we want ..

Any help is appreciated! I've also tried to solve it by integral transforms, but due to the stationarity (the problem-equation is homogenious), I fail as [itex] \rho [/itex] vanishes ..

 

[pasmith]

[itex]\rho = C[/itex] is a solution of the BVP. Have you any reason to expect others?

In general, what you want to look at in these sort of linear problems is [itex]u = \rho - C[/itex] rather than [itex]\rho[/itex] itself.

 

[mhsd91]

Indeed, you're absolutley right. My problem is that I'm working on a project where we're supposed to solve this eq. numerically, but by looking at it, I was really certain I could solve it analytically. It would be really nice to have such an analytical solution to compare with the numerical.

 

[eb90]

Maybe a Fourier transform could be a useful approach, as

[itex] \mathscr{F} \{{D\nabla^2 \rho(x,y) + \vec{V} \cdot \nabla \rho(x,y)}\}=-D(\mu^2+\nu^2) P(\mu,\nu)+\mathscr{F} \{ {\vec{V}} \}\ast i(\mu +\nu) P(\mu,\nu) [/itex],

or would this just give a condition for the vector field for the constant solution?

 

[mhsd91]

What you suggest is of course interesting. I fear however, I've already tried this: The Fourier (and Laplace) integral transform are indeed popular approaches for the non-stationary cases. Since I have stationarity, what you write there is equal to zero. Thus,

[itex]
0 = D (\mu^2 + \nu^2) P(\mu,\nu) + \mathscr{F}\{\vec{V}\} * i(\mu+\nu) P(\mu,\nu)
[/itex]

Obviously, we may divide this eq. by [itex] P(\mu,\nu) [/itex], and it is lost such that we cannot obtain [itex] \rho(x,y) [/itex] ..

Or am I doing something illegal here: will [itex] P(\mu,\nu) [/itex] be a part of the convolution? In that case, then a non constant solution of [itex] \rho(x,y) [/itex] may exist for non-constant [itex] \vec{V} [/itex]...?

 

[Chestermiller]

Indeed, you're absolutley right. My problem is that I'm working on a project where we're supposed to solve this eq. numerically, but by looking at it, I was really certain I could solve it analytically. It would be really nice to have such an analytical solution to compare with the numerical.

One thing you could do is to solve the problem for the case where the diffusion coefficient is zero. This could be done using the method of characteristics. This would give you a first order picture of what the solution for the density looks like. For small values of the diffusion coefficient, the solution would be a smoothed version of this.

Chet

 

[pasmith]

Indeed, you're absolutley right. My problem is that I'm working on a project where we're supposed to solve this eq. numerically, but by looking at it, I was really certain I could solve it analytically.

Indeed you can: set [itex]u = \rho - C[/itex]. Then [itex]u[/itex] satisfies [tex]D\nabla^2u + \mathbf{v} \cdot \nabla u = 0[/tex] with [itex]u(0,y) = u(x,0) = 0[/itex], [itex]\frac{\partial u}{\partial y} = 0[/itex] on [itex]y =1[/itex] and [itex]\frac{\partial u}{\partial x} = 0[/itex] on [itex]x = 1[/itex].

Separation of variables will work here: you'll have [itex]u = X(x)Y(y)[/itex] where [itex]X(0) = 0[/itex], [itex]X'(1) = 0[/itex], [itex]Y(0) = 0[/itex], and [itex]Y'(1) = 0[/itex].