Analytical approximations of the diffusive dispersion in fluid flows

Convection-diffusion equations occur in many physical and chemical problems and describe the transport of heat, matter, or other physical quantities due to a convective flow superimposed by the diffusion of the transported entity within the flow [1]. Only in exceptional cases is the convection-diffusion equation solvable by analytical means [2], but in most cases one has to take recourse to rather complicated numerical integration schemes. These numerical approaches include a wide array of sophisticated finite-element methods [3–12], lattice Boltzmann methods [13], particle methods [14], or solving the related Smoluchowski-Langevin stochastic differential equation [15], see also the overviews [16,17] and the excellent review [18]. In the present paper, we develop a path-integral–based approach to handling convection-diffusion equations and find a general approximate solution of the equation. We demonstrate our method on several examples, and consider in particular the diffusion within an advective shear flow profile over an interface. This problem is particularly important when dealing with the diffusion of molecules or particles close to the surface of microfluidic and nanofludic devices [19,20], or wetting and spreading phenomena [21].

The paper is organized as follows: In the next section, we present our general path-integral treatment of the diffusion-convection equation in two dimensions with an inhomogeneous flow along one direction (uni-directional flow). We obtain a general solution of the problem via a Fourier amplitude function that is expressed as a path integral involving the flow function. By using a second-order cumulant expansion, this path integral can be solved and yields a closed-form analytical approximation of the solution. Next, we apply our approximation to several examples and compare it with exact finite-element numerical solutions. Finally, we discuss our results and possible generalizations and extensions to more complex convection-diffusion problems.

In particular, we consider the convection-diffusion equation

$$\begin{equation} \frac{\partial c\left(\mathbf{r},t\right)}{\partial t} = D\Delta c\left( \mathbf{r}, t\right)-v\left(z\right) \partial_x c\left(\mathbf{r},t\right) \end{equation} \tag{ 1 }$$

where $c(\mathbf{r},t)$ is the local concentration of the transported and diffusing quantity at position r and time t, D is the diffusion coefficient, and v(z) is the z-dependent flow profile of a flow along the x-direction. For simplicity reasons, we restrict ourselves here to considering only two spatial dimensions, x and z, but the generalization to the third dimension will be straightforward. Also, we assume that the flow is directed only along the x-direction, but with an arbitrary flow profile v(z) as a function of the coordinate z. Green's function for this equation with the initial condition $G(\mathbf{r},\mathbf{r}_i,t=0)=\delta(\mathbf{r}-\mathbf{r}_i)$ can be written in path-integral form as [22]

$$\begin{equation} G(\mathbf{r}_f,\mathbf{r}_i,t) \propto \int\limits_{\left(\mathbf{r}_i, 0 \right)}^{\left(\mathbf{r}_f, t\right)}D\mathbf{r}(t^\prime)\exp\left[-S\left\{ \mathbf{r} (t^\prime)\right\} \right] \end{equation} \tag{ 2 }$$

where the functional S is defined by

$$\begin{equation} S = \int_0^t\text{d}t^\prime \frac{\left\{\dot{x}(t^{\prime})-v\left[ z\left(t^{\prime} \right)\right] \right\}^2 + \dot{z}^2(t^{\prime})}{4D}, \end{equation} \tag{ 3 }$$

where the path integration in eq. (2) runs over all possible trajectories starting at position $\mathbf{r}_i$ at time zero and ending at position $\mathbf{r}_f$ at time t. When considering a subset of all trajectories $\mathbf{r}(t')$ with one and the same $z(t')$, the path integration over $x(t')$ can be done analytically yielding the result

$$\begin{eqnarray} &\int\limits_{\left(x_i,0\right)}^{\left(x_f,t\right)} Dx(t^\prime) \exp\left[ - \int_0^t \text{d}t^{\prime} \frac{\left\{ \dot{x}(t^\prime) -v\left[z\left(t^\prime \right)\right] \right\}^2}{4D} \right] \rightarrow \nonumber\\ &\frac{1}{\sqrt{4 \pi D t}} \exp\left[-\frac{\left(x_f-x_i-\bar{x}\right)^2}{4 D t} \right], \end{eqnarray} \tag{ 4 }$$

where the quantity $\bar{x}$ is the flow transport along the x-direction measured along the specific trajectory $z(t')$ which is equal to

$$\begin{equation*} \int_0^t\text{d}t' v\left[z\left(t'\right)\right]. \end{equation*}$$

Thus, the final solution $G(\mathbf{r}_f, \mathbf{r}_i, t)$ is found to be

$$\begin{eqnarray} G(\mathbf{r}_f,\mathbf{r}_i,t)&=\int\text{d}\bar{x}\; w\left(\bar{x},z_f,z_i, t\right) \nonumber\\ &\quad \times\,\frac{1}{\sqrt{4\pi Dt}}\exp\left[-\frac{\left(x_f-x_i-\bar{x} \right)^2} {4Dt}\right],\quad \end{eqnarray} \tag{ 5 }$$

where the weight function $w\left(\bar{x},z_f,z_i,t\right)$ is given by the path integral

$$\begin{eqnarray} w\left(\bar{x},z_f,z_i,t\right)&\propto \int\limits_{\left(z_i,0\right)}^{\left( z_f, t\right)} Dz(t^\prime)\exp\left[ -\int_0^t\text{d}t^\prime \frac{\dot{z}^2(t^\prime)} {4D} \right]\nonumber\\ &\quad \times\,\delta\left\{\bar{x}-\int_0^t\text{d}t^\prime v\left[z(t^\prime)\right] \right\}. \end{eqnarray} \tag{ 6 }$$

Thus, in each plane $z=\text{const}$, the solution $G(\mathbf{r}_f,\mathbf{r}_i,t)$ is a superposition of Gaussian distributions with the distribution function $w(\bar{x},z_f,z_i,t)$ at center positions $\bar{x}$ along the x-axis. Now what remains is to find an analytical solution for the weight function $w(\bar{x},z_f, z_i,t)$. As we will see below, this can be done exactly for specific flow profiles v(z), but also in the general case of an arbitrary flow profile, the path integral of eq. (6) can be analytically approximated, which will be the topic of the next section.

To find an analytical approximation of eq. (6), we switch first to Fourier space by

$$\begin{equation} w\left(\bar{x},z_f,z_i,t\right)=\int\frac{\text{d}k}{2\pi}\tilde{w}\left(k,z_f,z_i,t \right) \exp(i k \bar{x}). \end{equation} \tag{ 7 }$$

Now, the Fourier amplitude $\tilde{w}$ is given by the path integral

$$\begin{eqnarray} &\tilde{w}\left(k,z_f,z_i,t \right) \propto \int\limits_{\left(z_i,0\right) }^{\left(z_f,t\right)} Dz(t^\prime) \nonumber\\ &\times\,\exp\left\{-\int_0^t\text{d}t^\prime \left[\frac{\dot{z}^2(t^\prime)}{4D} + i k v\left[z(t^\prime)\right]\right]\right\}. \end{eqnarray} \tag{ 8 }$$

To approximate this expression analytically, let us formally define the average over a functional $V[z(t')]$ by

$$\begin{eqnarray} &\left\langle V[z(t')] \right\rangle_{\left(z_i,0\right)}^{\left(z_f, t\right)}= \nonumber\\ &\frac{\int_{\left(z_i,0\right)}^{\left(z_f,t\right)} Dz(t^\prime) \exp \left\{-\int_0^t\text{d}t^\prime \frac{\dot{z}^2(t^\prime)}{4D}\right\} V[z(t')]} {\int_{\left(z_i,0\right)}^{\left(z_f,t\right)} Dz(t^\prime)\exp\left\{ -\int_0^t \text{d}t^\prime \frac{\dot{z}^2(t^\prime)}{4D} \right\}}. \end{eqnarray} \tag{ 9 }$$

The path integral in the denominator can be calculated analytically and is equal to

$$\begin{equation} g(z_f-z_i,t)=\frac{1}{\sqrt{4 \pi D t}}\exp\left[-\frac{\left(z_f-z_i\right)^2}{4 D t} \right]. \end{equation} \tag{ 10 }$$

Thus, the Fourier transform $\tilde{w}$ can be represented as an average over a functional in the form

$$\begin{eqnarray} \tilde{w}\left(k,z_f,z_i,t \right) &= g(z_f-z_i,t) \nonumber\\ &\quad \times\,\left\langle \exp\left\{-ik\int_0^t\text{d}t^\prime v\left[z(t^\prime) \right] \right\}\right\rangle_{\left(z_i,0 \right)}^{\left( z_f,t\right)}\qquad \end{eqnarray} \tag{ 11 }$$

The important next step is to apply a cumulant expansion to eq. (11) and to retain only the linear and the quadratic terms of the expansion:

$$\begin{eqnarray} &\left\langle\exp\left\{-ik\int_0^t\text{d}t^\prime v\left[z(t^\prime) \right]\right\} \right\rangle \approx \nonumber\\ &\,\exp\Biggl\{-ik\left\langle\bar{x}\right\rangle -\frac{k^2}{2} \left\langle \delta \bar{x}^2 \right\rangle\Biggr\}, \end{eqnarray} \tag{ 12 }$$

where the linear term is explicitly given by

$$\begin{eqnarray} &\left\langle\bar{x}\right\rangle_{z_f,z_i,t}=\left\langle\int_0^t\text{d}t^\prime v\left[ z(t^\prime)\right]\right\rangle_{\left(z_i, 0\right)}^{\left(z_f, t\right)} = \frac{1}{g(z_f-z_i,t)} \nonumber\\ &\times\,\int_0^t\text{d}t'\int\text{d}z'g(z_f-z',t-t')v\left(z'\right)g(z'-z_i,t') \qquad \end{eqnarray} \tag{ 13 }$$

and the quadratic term by

$$\begin{equation} \left\langle \delta \bar{x}^2 \right\rangle = \left\langle \bar{x}^2 \right\rangle - \left\langle \bar{x} \right\rangle^2, \end{equation} \tag{ 14 }$$

where the second moment is

$$\begin{eqnarray} \left\langle \bar{x}^2 \right\rangle_{z_f,z_i,t} &= \left\langle \int\limits_0^t \text{d}t^{\prime} \int\limits_0^t\text{d}t'' v\left[ z(t^\prime) \right] v\left[ z(t'') \right] \right\rangle_{z_i, 0}^{z_f, t} \nonumber\\ &= \frac{2}{g(z_f-z_i,t)}\int_0^t \text{d}t' \int\text{d}z' \nonumber\\ &\times g(z_f-z',t-t')v \left(z'\right) g(z'-z_i,t') \left\langle \bar{x} \right\rangle_{z',z_i,t'}.\nonumber\\ \end{eqnarray} \tag{ 15 }$$

With this cumulant approximation, eq. (12), one can back-transform $\tilde{w}$ into real space yielding an analytic expression for the weight function w:

$$\begin{eqnarray} &&w\left(\bar{x},z_f,z_i,t \right) \approx \frac{g(z_f-z_i,t)}{\sqrt{2 \pi \langle\delta\bar{x}^2\rangle}}\exp\left[-\frac{\left(\bar{x}-\langle\bar{x} \rangle\right)^2}{2 \langle\delta\bar{x}^2\rangle}\right].\quad \end{eqnarray} \tag{ 16 }$$

Then, the final approximation for Green's function of the convection-diffusion equation is found to be

$$\begin{eqnarray} &G(\mathbf{r}_f,\mathbf{r}_i,t)\approx \nonumber\\ &\frac{g\left(z_f-z_i,t\right)}{\sqrt{2\pi\left(2Dt+\langle\delta\bar{x}^2\rangle \right)}}\exp\left[-\frac{\left(x_f-x_i-\langle \bar{x} \rangle \right)^2}{4Dt+2 \langle\delta\bar{x}^2\rangle}\right],\quad \end{eqnarray} \tag{ 17 }$$

where $\langle \bar{x} \rangle$ and $\langle \delta \bar{x}^2 \rangle$ are themselves functions of z_i, z_f, and t as given by eqs. (13)–(15). Equation (17) is the core result of the present paper and approximates the solution of the convection-diffusion equation by a superposition of Gaussian distributions shifted along x by the mean value $\bar{x}$ which itself depends on the initial position z_i, and on the final position z_f where the solution is evaluated. Also, in addition to the spreading of the profile due to pure diffusion, 2Dt, there occurs a convection-induced spreading $\langle \delta \bar{x}^2 \rangle$ which also depends on z_i and z_f.

Let us first consider the special case of a linear flow profile $v(z) = v_1 z$ in free space. In that case, an exact solution is well known [23], and we will reproduce it here on the basis of eq. (17). For the linear profile, both the average $\langle\bar{x}\rangle_{z_f,z_i,t}$ and variance $\langle\delta\bar{x}^2\rangle_{z_f,z_i,t}$ can be calculated exactly:

$$\begin{equation} \langle\bar{x}\rangle_{z_f,z_i,t} = v_1 t \frac{z_i+z_f}{2}, \end{equation} \tag{ 18 }$$

$$\begin{equation} \langle\delta\bar{x}^2\rangle_{z_f,z_i,t}=\frac{Dv_1^2t^3}{6}. \end{equation} \tag{ 19 }$$

Moreover, for the linear profile, all higher cumulants are zero. Thus, in that case, the approximate solution given by eq. (17) is the exact solution. The left column of images in fig. 1 shows the temporal evolution of an initially Gaussian concentration profile with unity variance and center position $\{x,z\}=\{0,0\}$ for the time values $t=0,1,2,4$ and 8.

Zoom In Zoom Out Reset image size

Next, let us consider the parabolic flow profile $v(z) = v_2 z^2$ in free space. Again, both the average and variance can be calculated analytically yielding

$$\begin{eqnarray} \langle\bar{x}\rangle_{z_f,z_i,t} &= \frac{v_2 t (D t + z_i^2 + z_i z_f + z_f^2)}{3}, \end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray} \langle\delta\bar{x}^2\rangle_{z_f,z_i,t} &= \frac{2 D v_2^2 t^3 (2 D t + 4 z_f^2 +7 z_i z_f + 4 z_i^2)}{45}.\quad \end{eqnarray} \tag{ 21 }$$

However, now the higher-order cumulants do not vanish anymore, and eq. (17) is indeed only an approximation to the exact solution. The temporal evolution of this approximation is shown by the images in the middle column of fig. 1. To compare it with the exact solution, we integrated the full two-dimensional convection-diffusion problem with a finite-element numerical equation solver (COMSOL) on a grid stretching from $x=-80$ till x = 80 and $z=-20$ till z = 20 with a grid spacing of $\Delta x= \Delta z = 1/20$ and a time-step of 0.01. The full calculation between t = 0 and t = 8 lasted for 90 min on an Intel(R) Xeon(R) E5440 2,83 GHz quad core PC. For better visualization, we do not compare the full concentration distribution of the numerical calculation and the approximation, but the concentration profiles integrated along the z-axis. The result is shown in fig. 2. It should be noted that the concentration profiles integrated along the x-axis are identical because the integral of eq. (17) along x yields the exact result, independent of the particular function v(z). As can be seen, the approximation eq. (17) matches fairly well the dispersion of the concentration profile for small time values, and starts to deviate sensibly from the exact values only for large time values $t\gg1$.

Zoom In Zoom Out Reset image size

As a last example, we consider a linear flow profile over a planar surface with no-slip condition (flow velocity is zero on surface). This is the most relevant case for all convection-diffusion problems close to an interface. Instead of considering the convection-diffusion equation in a half-space z > 0 with flow profile $v(z)=v_1 z$, we will treat it in infinite space but assuming the cusp-like flow profile $v(z) = v_1 |z|$. When additionally imposing a mirror-symmetric initial distribution $c_0(x,z) = c_0(x,-z)$, the solution to this problem will be equivalent to the solution of the equation in half space over an impenetrable boundary. Now, there is no more an explicit analytical solution for the average $\langle\bar{x} \rangle_{z_f,z_i,t}$ and variance $\langle\delta\bar{x}^2\rangle_{z_f,z_i,t}$. Both quantities have to be computed numerically using eqs. (13)–(15). From eq. (13), we find the integral

$$\begin{equation} \left\langle\bar{x}\right\rangle_{z_f,z_i,t} =\int_0^t\text{d}t' \sqrt{\frac{4(t-t')t'}{t}} \left[\frac{e^{-\xi^2}}{\sqrt{\pi}}+\xi\text{erf} (\xi)\right], \end{equation} \tag{ 22 }$$

where we have introduced the abbreviation

$$\begin{equation} \xi = \frac{t' z_f +(t-t') z_i}{\sqrt{4 t (t-t') t'}}. \end{equation} \tag{ 23 }$$

This integral has to be evaluated numerically. Because the final result depends on three variables, namely z_f, z_i and time t, this could be a formidable task. However, the problem is much simplified when realizing that $\langle\bar{x} \rangle_{z_f,z_i,t}$ has the following time-scaling property:

$$\begin{equation} \left\langle\bar{x}\right\rangle_{z_f,z_i,t} = t^{3/2} \left\langle\bar{x} \right\rangle_{t^{-1/2} z_f,t^{-1/2} z_i,1}. \end{equation} \tag{ 24 }$$

Thus, one has to evaluate eq. (22) only for one specific value of time, t = 1, and can then find the function for all other values of time using eq. (24). Knowing $\langle\bar{x}\rangle_{z_f,z_i,t}$, one uses next the last equality in eq. (15) for numerically calculating the second moment $\langle\bar{x}^2\rangle_{z_f,z_i,t}$ and thus the variance $\langle \delta\bar{x}^2\rangle_{z_f,z_i,t}$. And similar to the average value $\langle \bar{x}\rangle_{z_f,z_i,t}$, one can use the time-scaling behavior of the variance,

$$\begin{equation} \left\langle \delta\bar{x}^2 \right\rangle_{z_f,z_i,t} = t^3 \left\langle \delta\bar{x}^2 \right\rangle_{t^{-1/2} z_f,t^{-1/2} z_i,1} \end{equation} \tag{ 25 }$$

for finding the function for arbitrary time values. The right column of images in fig. 1 shows the result for the temporal evolution of the Gaussian initial concentration profile for a linear cusp flow with $v_1=1$. The whole calculation took ca. 100 s on the computer. For comparison, we solved the full problem, as before, using the finite-element numerical equation solver COMSOL, where the full calculation again took ca. 90 min. Figure 3 shows a comparison of the temporal evolution of the z-integrated concentration profiles as obtained by full numerical integration and using the approximation (17). Again, the approximation yields fairly good results for small time values and starts to deviate visibly from the exact solution only for $t\gg1$.

Zoom In Zoom Out Reset image size

The core result of our paper is the approximate solution of the convection-diffusion equation (1) as given by eq. (17). The found approximation is very similar to the famous Taylor-Aris approximation for the convection-diffusion induced broadening of plug flow in capillaries [24,25]. There, the concentration profile along the flow direction is also approximated by a Gaussian distribution, for which the peak position and width are derived as functions of time. Furthermore, the physical meaning of the Gaussian approximation found in the present work is also similar to Taylor-Aris dispersion: The Gaussian approximation works very well as long as the diffusion transversal to the flow direction is fast enough for mixing the diffusing molecules perpendicular to that direction before the flow gradient disrupts the profile. To quantify this, one introduces the Péclet number, $\text{Pe} = \bar{v} L/D$, where $\bar{v}$ is the mean flow velocity, L the transversal size of the capillary, and D the diffusion coefficient, which describes the ratio of diffusive to convective transport. As long as the Péclet number is small, i.e. as long as the diffusive transport perpendicular to the flow direction dominates over the convective transport, the Taylor-Aris approximation works well. In our examples considered in the present letter, there is no intrinsic lateral size L, because the considered systems are unbounded. However, we consider the transport of an initial Gaussian distribution with distribution width $\sigma_0$. The width of this distribution perpendicular to the flow direction evolves as $\sigma_\perp(t) = \sqrt{\sigma_0^2+2 D t}$, and we can take $\sigma_\perp(t)$ as the momentous characteristic transversal size of the problem. Then, the average flow velocity over this length scale is proportional to $v_1 \sigma_\perp(t)$ for the linear cusp velocity profile, and to $v_2 \sigma_\perp(t)^2$ for the parabolic velocity profile. Thus, one can define the time-dependent Péclet number $\text{Pe}_1(t) = v_1 \sigma_\perp(t)^2/D$ and $\text{Pe}_2(t) = v_2 \sigma_\perp(t)^3/D$ for the linear and for the parabolic flow profile problems, respectively. As long as these numbers are less or on the order of unity, the Gaussian approximation will work well. But this simple analysis also shows that the fidelity of the approximation deteriorates with time approximately proportional to t for the case of the linear cusp profile, and proportional to $t^{3/2}$ for the case of the parabolic profile. The main reason is that we considered here only unbounded problems. For a bounded system, where diffusion will be limited in space perpendicular to the flow direction, the approximation will work over longer time-scales, similar to the Taylor-Aris case with low Péclet numbers.

This leads to another important discussion: The approach presented in this letter can be easily generalized to any problem where the convection-diffusion problem can be formulated as

$$\begin{equation} \frac{\partial c}{\partial t} = D \left(\Delta_\perp + \frac{\partial^2}{\partial x^2}\right)c+v(\mathbf{r}_\perp)\frac{\partial c}{\partial x}, \end{equation} \tag{ 26 }$$

where $\{x,\mathbf{r}_\perp\}$ form an orthogonal coordinate system and $\Delta_{\perp}$ being the part of the Laplace operator not containing any derivative with respect to x. Then the function g in eqs. (13), (15) is the fundamental solution of the transverse diffusion problem

$$\begin{equation} \frac{\partial g}{\partial t} = D \Delta_\perp g \end{equation} \tag{ 27 }$$

and the integration in these equations extends over all transverse coordinates $\mathbf{r}_\perp$. Importantly, this concept is also applicable to transversally bounded problems such as shear flow between two parallel plates moving laterally with respect to each other.

Finally, we want to emphasize that the new path integral representation of the solution as given by eqs. (5), (7) and (8) can be used also for finding other exact solutions of convection-diffusion problems, in particular for all cases where the path integral (8) can be evaluated analytically, or for devising efficient numerical integration schemes.

Financial support by the Deutsche Forschungsgemeinschaft is gratefully acknowledged (SFB 937, project A5). NK thanks the Niedersächsische Ministerium für Wissenschaft und Kultur for an MWK stipend.