# Diluted-dispersed mass transfer within an AWE

N.Valle  Process & Energy Department, Delft University of Technology,
Leegwaterstraat 39, 2628 CB Delft, The Netherlands

###### Abstract

The goal of this document is describe the multiphase transfer processes
describing the bubble dynamics of a water electrolyzer. The motivation is to
describe the dilute-dispersed mass transfer within and Alkaline Water
Electrolyzer. Special emphasis is put on the mathematical formulation. The
presentation starts by posing the governing equations and their dimensionless
counterpart. By filtering the equations, the two-fluid model is presented
along with the need to sub-scale and wall models. To the later aim, boundary
layer equations are introduced. By reviewing self-similiarity transformations,
the analysis of Blasius, Ostrach and Sparrow is reviewed for Prandtl’s
boundary layer equations; along with that of Leveque.

##  1 Governing equations

Conservation laws

We begin with the formulation of mass, momentum and species conservation:

$\displaystyle\frac{\partial\rho}{\partial t}+\nabla\cdot\rho\vec{u}$
$\displaystyle=0$ (1) $\displaystyle\frac{\partial\rho\vec{u}}{\partial
t}+\nabla\cdot\rho\vec{u}\otimes\vec{u}$ $\displaystyle=-\nabla
p+\nabla\cdot\tau+\rho\vec{g}$ (2) $\displaystyle\frac{\partial c}{\partial
t}+\nabla\cdot c\vec{u}$ $\displaystyle=D\nabla\cdot f+kc$ (3)

Boundary conditions

Boundaries consist of the channel walls along with the input and output
regions of the channel.

$\displaystyle p\rvert_{i}$ $\displaystyle=p_{0}$ $\displaystyle p\rvert_{i}$
$\displaystyle=0$ $\displaystyle\vec{u}\rvert_{o}\hat{t}_{o}$
$\displaystyle=0$ $\displaystyle\vec{u}\rvert_{w}$ $\displaystyle=0$ (4)
$\displaystyle c\rvert_{i}$ $\displaystyle=c_{s}$ $\displaystyle D\nabla
c\rvert_{o}\cdot\hat{t}_{o}$ $\displaystyle=0$ $\displaystyle D\nabla
c\rvert_{w}\cdot\hat{n}_{w}$ $\displaystyle=kc$ (5)

Local instant formulation

Now, introducing the marker function $\phi$, we proceed to the splitting of
the continuum equations between continuum (i.e., liquid) and dispersed (i.e.,
gas) phases.

$\displaystyle\frac{\partial\phi_{k}\rho}{\partial
t}+\nabla\cdot\phi_{k}\rho\vec{u}$ $\displaystyle=m_{k}$ (6)
$\displaystyle\frac{\partial\phi_{k}\rho\vec{u}}{\partial
t}+\nabla\cdot\phi_{k}\rho\vec{u}\otimes\vec{u}$ $\displaystyle=-\nabla\phi
p+\mu\nabla\cdot\phi_{k}S+\phi_{k}c\vec{g}+M_{k}$ (7)
$\displaystyle\frac{\partial\phi_{k}c}{\partial
t}+\nabla\cdot\phi_{k}c\vec{u}$
$\displaystyle=D\nabla^{2}\phi_{k}c+k\phi_{k}c$ (8)

where $m_{k}$ corresponds to the mass transfer between phases, while $M_{k}$
corresponds to the momentum transfer.

$\displaystyle\dot{m}_{k}$
$\displaystyle=\rho_{k}(\vec{u}_{\Gamma}-\vec{u})\cdot\hat{n}_{\Gamma}\delta_{\Gamma}$
(9) $\displaystyle M_{k}$
$\displaystyle=\dot{m}_{k}\cdot\vec{u}+\sigma\hat{n}\delta_{\Gamma}$ (10)

Incompressible

$\vec{u}=0$ (11)

Newtonian

We will assume that the fluid is Newtonian, and so that the relationship
between strain ($S$) and stress ($\tau$) is linear:

$\tau=\mu S$ (12)

Boussinesq While both density and viscosity are functions of the local
supersaturation, we will adopt here the Boussinesq hypothesis. This is
reasonable as far as the concentrations are close to saturation. Consequently,
we will consider them constant everywhere except for the formulation of the
buoyancy term, where we approximate it with a 2nd order Taylor expansion:

$\rho(c)\approx\rho_{s}+\frac{\partial\rho}{\partial c}(c-c_{s})$ (13)

Nonetheless, since the dispersed phase is composed of very small bubbles,
whose concentration gradients within the bubble itself are expected to be very
small, buoyancy effects inside the bubbles have little room to develop and so
are neglected. Actually, we consider concentration to be constant within the
dispersed phase and equal to $c_{s}$. Note, however, that density differences
across the bubble-liquid interface do induce strong buoyancy effects on the
bubbles relative to the liquid phase, which are effectively included.
Consequently, we will only introduce the Boussinesq hypothesis for the liquid
phase (i.e., $k=l$).

Final form

By referencing the buoyancy terms to $\rho_{s}$, we finally obtain:

$\displaystyle\frac{\partial\phi_{k}}{\partial t}+\nabla\cdot\phi_{k}\vec{u}$
$\displaystyle=\dot{m}_{k}$ (14)
$\displaystyle\frac{\partial\phi_{l}\rho\vec{u}}{\partial
t}+\nabla\cdot\phi_{l}\rho\vec{u}\otimes\vec{u}$
$\displaystyle=-\nabla\phi_{l}p+\nabla\cdot\phi_{l}\mu
S+\phi_{l}\frac{\partial\rho}{\partial c}(c-c_{s})\vec{g}+M_{l}$ (15)
$\displaystyle\frac{\partial\phi_{g}\rho\vec{u}}{\partial
t}+\nabla\cdot\phi_{g}\rho\vec{u}\otimes\vec{u}$
$\displaystyle=-\nabla\phi_{g}p+\nabla\cdot\phi_{g}\mu
S+\phi_{g}(\rho-\rho_{l})\vec{g}+M_{g}$ (16)
$\displaystyle\frac{\partial\phi_{l}c}{\partial t}+\nabla\cdot\phi_{l}cu^{i}$
$\displaystyle=D\nabla^{2}\phi_{l}c+k\phi_{l}c-\dot{s}_{l}$ (17)
$\displaystyle\phi_{g}c$ $\displaystyle=c_{s}$ (18)

which consitutes the local instantious formulation of the governing equations.
We can see that liquid and gas phases exhibit a different nature of the
buoyancy terms (concentration driven in the liquid, vs. density driven in the
gas). We also see that species transport is only relevant within the liquid
phase. Note the difference between the effect between species (and mass)
transfer in both dispersed and continuos phases: for the gas phase, whose
concentration is constant, $\dot{s}_{g}$ contributes only to the expansion of
the phase via the conservation of mass; whereas for the liquid phase, whose
concentration is variable, $\dot{s}_{l}$ both contributes to the shrinking of
the gas phase and its dilution.

###  1.1 Modeling of closure terms

Mass transfer coefficient

While the physical mechanism of mass transfer is described by the diffusion of
dissolved gas from the liquid into the bubble, it is most of the time
convenient to model such a flow by a mass transfer coefficient $k_{b}$

$D\nabla c\rvert_{\Gamma}=k_{b}(c-c_{s})\rvert_{\Gamma}$ (19)

such that we obtain the flux of species as:

$\dot{s}_{k}=\pm k_{b}(c-c_{s})\delta_{\Gamma}$ (20)

which turns into a mass flux by using the molar density $N$

$\dot{m}_{k}=\pm Nk_{b}(c-c_{s})\delta_{\Gamma}$ (21)

##  2 Dimensionless formulation

We now introduce the general dimensionless variables for a buoyancy driven
flow:

$\displaystyle x$ $\displaystyle=LX$ $\displaystyle t$
$\displaystyle=t_{ref}\tau=\frac{\rho L^{2}}{\mu}\tau$ $\displaystyle u$
$\displaystyle=\frac{L}{t_{ref}}U=\frac{\mu}{\rho L}U$ $\displaystyle p$
$\displaystyle=\rho_{0}u_{ref}^{2}P=\frac{\mu^{2}}{\rho L^{2}}P$
$\displaystyle c$ $\displaystyle=c_{s}\left(\zeta+1\right)$
$\displaystyle\rho$ $\displaystyle=\rho_{0}\rho^{*}$ $\displaystyle\mu$
$\displaystyle=\mu_{0}\mu^{*}$ $\displaystyle D$ $\displaystyle=D_{0}D^{*}$

where we have used phase $0$ as the reference values of the physical
properties, and introduced $\rho^{*}$, $\mu^{*}$ and $D^{*}$ as dimensionless
properties.

With this formulation, we can restate mass, species and momentum conservation
in dimensionless form as:

$\displaystyle\frac{\partial\phi_{k}}{\partial t}+\nabla\cdot\phi_{k}\vec{u}$
$\displaystyle=\pm\boxed{\frac{Nc_{s}}{\rho_{k}}}~{}\boxed{\frac{k_{b}L\rho_{0}}{\mu}}\zeta\delta_{\Gamma}$
(22) $\displaystyle\frac{\partial\phi U^{i}}{\partial t}+\nabla\cdot\phi
U^{i}U^{j}$ $\displaystyle=-\nabla\phi P+\nabla\cdot\phi
S+\boxed{\frac{gL^{3}c_{s}\rho_{s}}{\mu_{s}}\frac{\partial\rho}{\partial
c}}\zeta\vec{g}$ (23) $\displaystyle\frac{\partial\phi\rho^{*}U^{i}}{\partial
t}+\nabla\cdot\phi\rho^{*}U^{i}U^{j}$ $\displaystyle=-\nabla\phi
P+\nabla\cdot\phi\mu^{*}S+\boxed{\frac{\rho_{l}(\rho_{g}-\rho_{l})gL^{3}}{\mu_{0}^{2}}}\rho^{*}\zeta\vec{g}$
(24) $\displaystyle\frac{\partial\phi\zeta}{\partial t}+\nabla\cdot\zeta
U^{i}$
$\displaystyle=\boxed{\frac{D\rho_{s}}{\mu_{s}}}D^{*}\nabla^{2}\zeta+\boxed{\frac{kL^{2}\rho_{s}}{\mu_{s}}}(\zeta+1)+\boxed{\frac{k_{b}L\rho}{\mu}}\zeta$
(25)

and also the boundary conditions

$\displaystyle\frac{\mu}{\rho L}U\rvert_{wall}$ $\displaystyle=0$ (26)
$\displaystyle\frac{Dc_{s}}{L}\nabla\zeta\rvert_{wall}\cdot\hat{n}_{wall}$
$\displaystyle=kc_{s}(\zeta+1)$ (27)

Which can be rewritten most compactly by identifying the characteristic
dimensionless numbers as follows:

$\displaystyle\frac{\partial\phi_{l}}{\partial t}+\nabla\cdot\phi_{l}U$
$\displaystyle=-\frac{\rho_{g}}{\rho_{l}}\frac{Sh_{b}}{Sc}\zeta\delta_{\Gamma}$
(28) $\displaystyle\frac{\partial\phi_{g}}{\partial t}+\nabla\cdot\phi_{g}U$
$\displaystyle=+\frac{Sh_{b}}{Sc}\zeta\delta_{\Gamma}$ (29)
$\displaystyle\frac{\partial\phi U}{\partial t}+\nabla\cdot\phi U\otimes U$
$\displaystyle=-\nabla\phi P+\nabla\cdot\phi S+Gr\zeta\hat{g}+M_{l}$ (30)
$\displaystyle\frac{\partial\phi\rho^{*}U}{\partial
t}+\nabla\cdot\phi\rho^{*}U\otimes U$ $\displaystyle=-\nabla\phi
P+\nabla\cdot\phi\mu^{*}S+\rho^{*}Ar\hat{g}+M_{g}$ (31)
$\displaystyle\frac{\partial\phi\zeta}{\partial t}+\nabla\cdot\zeta U$
$\displaystyle=\frac{1}{Sc}\nabla^{2}\zeta+\frac{Da_{\Pi}}{Sc}(\zeta+1)+\frac{Sh_{b}}{Sc}\zeta$
(32) $\displaystyle\phi_{g}\zeta$ $\displaystyle=1$ (33)

subject to

$\displaystyle U\rvert_{wall}$ $\displaystyle=0$ (34)
$\displaystyle\nabla\zeta\rvert_{wall}\cdot\hat{n}_{wall}$
$\displaystyle=Da(\zeta+1)$ (35)

where $Gr$ is the Grashoff number, $Ar$ is the Archimedes number, $Sc$ is the
Schmidt number, $Da$ is the (first) Damkhöler number, and $Da_{\Pi}$ is the
second Damkhöler number.

This body of equations rules the overall mass transfer phenomena occurring in
a dissolved/dispersed system. Let’s explore the role of several terms in
further detail: The Grashoff number ($Gr$) represents the buoyancy effects of
the liquid phase due to the density gradients induced by the concentration
gradients. However, this buoyancy effect is rapidly exceed by bubble-induced
buoyancy. This is produced by the density difference between the bubble and
its surrounding liquid and fueled by gravity. The intensity of this effect is
represented by the Archimedes ($Ar$) number.

##  3 Filtering

The solution of the aforementioned equations is know to exhibit very fine
spatio-temporal features. These are mainly due to two main physical process:
(i) the chaotic (turbulent) nature of the solution, but also due to  (ii)
small-scale flow dynamics  However, in most engineering applications we are
interested in the large scale features. Mathematically, we can formulate these
large scale features as a filtering of the solution. We are then interested in
obtaining filtered form of the governing equations, which (ideally) is
expressed in terms of the filtered variables only.

To do so, we introduce the filtering operator

$\langle\cdot\rangle=\frac{1}{\Omega}\int_{\Omega}\cdot d\Omega$ (36)

where $\Omega$ could be a characteristic time, space or thermodynamic state,
giving rise to time- space- or ensemble-averaged filters. After filtering, one
can perform a Reynolds decomposition of any variable as an average value plus
an oscillation from that average

$u=\overline{u}+u^{\prime}$ (37)

subject to $\langle u^{\prime}\rangle=0$.

By filtering the governing equations 38-43, we obtain:

$\displaystyle\frac{\partial\epsilon_{l}}{\partial
t}+\nabla\cdot\epsilon_{l}\tilde{U}$
$\displaystyle=-\frac{\rho_{g}}{\rho_{l}}\frac{Sh_{b}}{Sc}\epsilon_{l}\overline{\zeta}\delta_{\Gamma}-\boxed{\nabla\cdot\overline{\phi_{l}^{\prime}u^{\prime}}}$
(38) $\displaystyle\frac{\partial\epsilon_{g}}{\partial
t}+\nabla\cdot\epsilon_{g}\tilde{U}$
$\displaystyle=+\frac{Sh_{b}}{Sc}\epsilon_{g}\tilde{\zeta}\delta_{\Gamma}-\boxed{\nabla\cdot\overline{\phi_{g}^{\prime}u^{\prime}}}$
(39) $\displaystyle\frac{\partial\epsilon_{l}\tilde{U}}{\partial
t}+\nabla\cdot\epsilon_{l}\tilde{U}\otimes\tilde{U}$
$\displaystyle=-\nabla\epsilon_{l}\overline{P}+\nabla\cdot\epsilon_{l}\tilde{S}+Gr\overline{\zeta}\hat{g}+\overline{M_{l}}-\boxed{\nabla\cdot\epsilon\overline{u^{\prime}\otimes
u^{\prime}}}$ (40)
$\displaystyle\rho^{*}\frac{\partial\epsilon_{g}\tilde{U}}{\partial
t}+\nabla\cdot\epsilon_{g}\tilde{U}\otimes\tilde{U}$
$\displaystyle=-\nabla\epsilon_{g}\overline{P}+\mu^{*}\nabla\cdot\epsilon_{g}\tilde{S}+\rho^{*}Ar\hat{g}+\overline{M_{g}}-\boxed{\nabla\cdot\epsilon\overline{u^{\prime}\otimes
u^{\prime}}}$ (41)
$\displaystyle\frac{\partial\epsilon_{l}\overline{\zeta}}{\partial
t}+\nabla\cdot\epsilon_{l}\overline{\zeta}\tilde{U}$
$\displaystyle=\frac{1}{Sc}\nabla^{2}\epsilon_{l}\overline{\zeta}+\frac{Da_{\Pi}}{Sc}\epsilon_{l}(\overline{\zeta}+1)+\frac{Sh_{b}}{Sc}\epsilon_{l}\overline{\zeta}-\boxed{\nabla\cdot\epsilon\overline{\zeta^{\prime}u^{\prime}}}$
(42) $\displaystyle\epsilon_{g}\overline{\zeta}$ $\displaystyle=1$ (43)

where we have introduced new terms (boxed) arising from the oscillatory
behavior of the solution. These terms require a sub-grid scale model in order
to close the relations.

Let’s analyze them in detail:

  * $\overline{u^{\prime}\otimes u^{\prime}}$

this terms is the Reynolds stress tensor and describes the additional stress
experienced by the flow due to the turbulent oscillations. The most common way
to model these terms are the eddy viscosity models, which assimilate them with
an additional diffusion term.

  * $\overline{\zeta^{\prime}u^{\prime}}$

this terms is analogous to the turbulent dispersion force. However, it has two
contributions: the purely turbulent dispersion, and the effect of the
micromixing produced in the vicinity of the wall. Assuming the Reynolds
analogy, we can assume that the turbulent dispersion is proportional to the
eddy viscosity, and so can be seen as the modification of the local Schmidt
number. Regarding the micromixing effects, we need to come up with an
additional closure term to model them.

###  3.1 Modeling of sub-scale terms

Mass transfer

We will assess here the contribution of sub-scale terms to the mass transfer
equation. In particular, we will look at the filtered term of equation (43)
$\overline{\zeta^{\prime}\otimes u^{\prime}}$. This term includes two
different contributions which are very different in its nature.

The first contribution is micromixing, which is due to a perturbation within
the mass transfer boundary layer. When the solution is not locally isotropic
within a computational domain (e.g., near a wall) filtering will result smear
the solution, introducing a non-physical damping of the solution. In this
cases, we need to advance the resolution of the computation by so-called wall
models.

The second contribution is turbulent dispersion due to the turbulent
oscillations of the system. On the other hand, even when the solution is
oscillatory within the computational domain, since non–linear terms (e.g.,
convection) do not commute with filtering, result into additional terms due to
the propagation of the oscillations. Such additional terms need to be modeled
by what we call sub-grid scale models. This term can be modeled analogously to
the Reynolds stress term reported before by means of and eddy viscosity model.

Wall model Since the typical size of mass transfer boundary layers is much
smaller than the hydrodynamic ones, since they scale as:

$\frac{\delta}{\delta_{c}}=Sc^{-1/3}$ (44)

for typical values of Schmidt numbers $Sc>>1$ we obtain $\delta_{c}<<\delta$.
Since the resolution required to resolve $\delta$, the hydrodynamic boundary
layer, is already very high, we can model the mass transfer boundary by means
of a Sherwood-like expression as:

$\nabla\zeta\rvert_{w}=Sh_{w}c$ (45)

which allows to work with a much coarser mesh.

$0$$0.2$$0.4$$0.6$$0.8$$1$$0$$0.5$$1$$1.5$$\delta_{c}$$\delta$$c_{w}$$c_{i}$$\Delta
x$$c_{i}-c_{w}$$x$ Figure 1: The concentration profile around the eletrode
wall exhibits a really thin boundary layer ($\delta_{c}$), which is much
thinner than the Blasius hydrodynamic boundary layer ($\delta$), and the
characteristic filter length ($\Delta x$). Obtaining the diffusive flow by a
classical finite difference scheme will result into an artificially low
diffusion flux, since $\Delta c/\Delta x<<\Delta c/\delta_{c}$, so we
typically tune the value of $D$ by means of a wall model.

While a number of relations for $Sh$ exist for vertical planes, the presence
bubbles detaching from the electrode surface intensify the mass transfer of
species towards the bulk. Jansen et al. (1989) reported on the effect of gas
bubbles, distinguishing between two main mechanisms of bubble detachment:  (i)
bubble sliding, detachment and rising; and  (ii) bubble sliding, coalescing
and jumping.  The former corresponds with bubbles that typically do not
coalesce, like it has been reported for $H_{2}$, which slide due to gravity
over the electrode wall before departing into the bulk. The later corresponds
with bubbles that do nucleate, like it has been reported for $O_{2}$, which
tend to coalesce violently on the electrode surface and violently jump into
the bulk of the channel.

The data reported there can expressed as the superposition of the average $Sh$
number at a bubble-free wall, plus an enhancement factor produced by the
bubbles. This model has obtained an outstanding match with experimental
results, as can be seen in Janssen1989.

We obtain, for sliding and rising bubbles:

$Sh_{w,b}=Sh_{w}\left(1+v_{g\to
b}\left(\frac{6\alpha_{1}^{3}s}{4\pi^{2}D}\right)\right)+\frac{3\alpha_{2}^{2}R}{4\pi
DSh_{w}}$ (46)

Sub-scale model

##  4 Boundary layer

Boundary layers occur in the vicinity of the wall, where viscosity is dominat.
This allows to take several simplifications on the governing equations, which
enables some analytical incursion in the flow solutions.

$0$$0.2$$0.4$$0.6$$0.8$$1$$0$$0.5$$1$$1.5$$\delta_{c}$$\delta$$x$ Figure 2:
The concentration profile around the eletrode wall exhibits a really thin
boundary layer ($\delta_{c}$), which is much thinner than the Blasius
hydrodynamic boundary layer ($\delta$).

###  4.1 Prandtl boundary layer equations

Ludwig’s Prandtl model for the boundary layer equations assumes a 2D,
incompressible, stationary flow in the vicinity of a plate. This results into
the following set of equations:

$\displaystyle\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}$
$\displaystyle=0$ (47) $\displaystyle u\frac{\partial u}{\partial
x}+v\frac{\partial u}{\partial y}$
$\displaystyle=-\frac{1}{\rho}\frac{\partial p}{\partial
x}+\nu\frac{\partial^{2}u}{\partial y^{2}}$ (48) $\displaystyle\frac{\partial
p}{\partial y}$ $\displaystyle=0$ (49)

assume that the variations along the stream-wise direction are much smaller
than the variation in the wall-normal one.

Note that equations (47)-(49) involve three unknowns ($u$, $v$ and $p$) and
three equations.

To simplify the equation a little bit further, we may assume that pressure
variation in the streamwise direction is either zero (as for the case of
boundary plates) or just hydrostatic (i.e., $p=\rho gx$).

However, since the flow is incompressible, we can express velocities in terms
of a single streamfunction, $\Psi$.

$\displaystyle u$ $\displaystyle=\frac{\partial\Psi}{\partial y}$
$\displaystyle v$ $\displaystyle=-\frac{\partial\Psi}{\partial x}$ (50)

This results in the following, single-variable PDE:

$\displaystyle\frac{\partial\Psi}{\partial y}\frac{\partial^{2}\Psi}{\partial
x\partial y}-\frac{\partial\Psi}{\partial x}\frac{\partial^{2}\Psi}{\partial
y\partial y}$ $\displaystyle=\nu\frac{\partial^{3}\Psi}{\partial y^{2}y}+b$
(51)

with the following boundary conditions:

$\displaystyle\left.u\right|_{y=0}$
$\displaystyle=\left.\frac{\partial\Psi}{\partial y}\right|_{y=0}=0$
$\displaystyle\left.v\right|_{y=0}$
$\displaystyle=-\left.\frac{\partial\Psi}{\partial x}\right|_{y=0}=0$
$\displaystyle\left.u\right|_{x\to\infty}$
$\displaystyle=\left.\frac{\partial\Psi}{\partial x}\right|_{x\to\infty}=0$
(52)

###  4.2 Transport of scalars

The mass transport equation of a scalar within the boundary layer is:

$u\frac{\partial\theta}{\partial x}+v\frac{\partial\theta}{\partial
y}=D\frac{\partial^{2}\theta}{\partial y^{2}}$ (53)

which can be expressed in terms of the streamfunction as:

$\frac{\partial\Psi}{\partial y}\frac{\partial\theta}{\partial
x}-\frac{\partial\Psi}{\partial x}\frac{\partial\theta}{\partial
y}=D\frac{\partial^{2}\theta}{\partial y^{2}}$ (54)

with the following boundary conditions

$\displaystyle\left.\Delta c\right|_{y=0}$ $\displaystyle=c_{wall}-c_{\infty}$
$\displaystyle\left.\Delta c\right|_{y\to\infty}$ $\displaystyle=0$ (55)
$\displaystyle\left.D\Delta c\right|_{y=0}$ $\displaystyle=q$
$\displaystyle\left.\Delta c\right|_{y\to\infty}$ $\displaystyle=0$ (56)

for the Dirichlet or the Neumann case, respectively.

##  5 Self-similarity solutions

Equation (51) involves two variables ($x$ and $y$), which complicates its
solution by analytical or numerical methods.

Self-similarity solutions are solutions to PDEs in terms of a self-similarity
variable, $\eta$ which in general it is a function of both $x$ and $y$. In
this way, we can say that the solutions of the system at any point in space
(i.e., in terms of $x$ and $y$) are “similar” to the solution in terms of
$\eta$.

This tool allows to reduce the number of variables, provided that a self-
similarity variable can be found or guessed.

Since most of the times we do not know the exact form of the self-similarity
function _a priori_ , we attempt generic functions for the self-similarity
variables and functions, substitute them in the formulation of the original
PDE and finally the generic function to make sure that the resulting PDE is a
function of the self-similarity variables solely. We typically do so by
assuming some sort of polynomial function on $x$ and $y$. The procedure is the
following:

First, we do assume a shape for the self-similarity variable $\eta$ and the
self-similarity solution $\Psi$. These functions can, in principle, be any
analytic function. As fluid dynamicists, this is where our understanding of
the physics of the problem, funneled trough scaling arguments, may lead to
fine educated guesses on the shape of the equations.

On the one hand, we assume that the solution is self-similar in the streamwise
direction. This can be hinted from the shape of the equations, since it
includes a purely convective term (which has a purely transport character) in
the streamwise direction, while diffusion (which has a spreading character) in
the wall-normal direction. Since we are looking at boundary layer flows, we
render $\eta=y/\delta$ to represent the ratio of the wall-normal direction $y$
to the thickness of the boundary layer $\delta$, being $\delta$ some function
of the streamwise coordinate $x$. We further hint, that $\delta$ is
proportional to some power of $x$. As such, our candidate $\eta$ looks like:

$\eta=Ayx^{-s}$ (57)

The velocity field is self-similiar in the streamwise direction. However,
since velocity is the derivative of the streamfunction, $\Psi$, this implies
that $\Psi$ is a combination of a self-similar function $f(\eta)$ and some
power of $x$. This renders our candidate $f$ as:

$f(\eta)=B\Psi x^{-t}$ (58)

The actual values of $s$ and $t$ is what gives _the scalings_ of the solution.

While some candidates may exist for $s$ and $t$ from experimental fittings, we
cannot say much on its value _a priori_. As such, we will then proceed to
obtain a self-similarity solution for equation (51).

We start by describing the following intermediate variables:

$\displaystyle u$ $\displaystyle=\frac{\partial\Psi}{\partial
y}=\frac{\partial\eta}{\partial y}\frac{\partial\Psi}{\partial\eta}$
$\displaystyle=\frac{A}{B}x^{t-s}f^{\prime}$ $\displaystyle-v$
$\displaystyle=\frac{\partial\Psi}{\partial x}=\frac{\partial\eta}{\partial
x}\frac{\partial\Psi}{\partial\eta}$
$\displaystyle=\frac{1}{B}x^{t-1}\left(tf-s\eta f^{\prime}\right)$
$\displaystyle\frac{\partial u}{\partial x}$
$\displaystyle=\frac{\partial^{2}\Psi}{\partial x\partial y}$
$\displaystyle=\frac{A}{B}x^{t-s-1}\bigl{(}(t-s)f^{\prime}-s\eta
f^{\prime\prime}\bigr{)}$ $\displaystyle\frac{\partial u}{\partial y}$
$\displaystyle=\frac{\partial^{2}\Psi}{\partial y^{2}}$
$\displaystyle=\frac{A^{2}}{B}x^{t-2s}f^{\prime\prime}$
$\displaystyle\frac{\partial^{2}u}{\partial y^{2}}$
$\displaystyle=\frac{\partial^{3}\Psi}{\partial y^{3}}$
$\displaystyle=\frac{A^{3}}{B}x^{t-3s}f^{\prime\prime\prime}$

where we have adopted the prime notation $(^{\prime})$ to denote
differentiation with respect to $\eta$ for the sake of compactness. Since $f$
is a function of $\eta$ only, and all differentiation is taken with respect to
$\eta$ solely, there is no risk of confusion.

This produces:

$f^{\prime\prime\prime}+\frac{1}{\nu
AB}x^{t+s-1}\left(tff^{\prime}-(t-s)f^{\prime 2}\right)=b\frac{B}{\nu
A^{3}}x^{3s-t}$ (59)

with the following boundary conditions:

$\displaystyle\left.u\right|_{y=0}$
$\displaystyle=0=\left.\frac{A}{B}x^{t-s}f^{\prime}\right|_{\eta=0}$
$\displaystyle\left.f^{\prime}\right|_{0}$ $\displaystyle=0$ (60)
$\displaystyle\left.v\right|_{y=0}$
$\displaystyle=0=-\left.\frac{1}{B}x^{t-1}\left(tf-s\eta
f^{\prime}\right)\right|_{\eta=0}$ $\displaystyle\left.f\right|_{0}$
$\displaystyle=0$ (61) $\displaystyle\left.u\right|_{x\to\infty}$
$\displaystyle=0=\left.\frac{A}{B}x^{t-s}f^{\prime}\right|_{x\to\infty}$
$\displaystyle\left.f^{\prime}\right|_{\infty}$ $\displaystyle=0$ (62)

Introducing the scalar self-similarity variable:

$\theta(\eta)=C\Delta cx^{-u}$ (63)

introducing the following intermediate variables

$\displaystyle\frac{\partial c}{\partial x}$
$\displaystyle=\frac{1}{C}x^{u-1}\bigl{(}u\theta^{\prime}-s\eta\theta^{\prime}\bigr{)}$
$\displaystyle\frac{\partial c}{\partial y}$
$\displaystyle=\frac{A}{C}x^{u-s}\theta^{\prime}$
$\displaystyle\frac{\partial^{2}c}{\partial y^{2}}$
$\displaystyle=\frac{A^{2}}{C}x^{u-2s}\theta^{\prime\prime}$

and replacing them in the streamfunction formulation of the mass transport
equation (53), we obtain:

$\theta^{\prime\prime}+\frac{1}{DAB}x^{t+s-1}\left(tf\theta^{\prime}-uf^{\prime}\theta\right)=0$
(64)

with the corresponding boundary conditions:

$\displaystyle\left.\theta\right|_{0}$ $\displaystyle=1$
$\displaystyle\left.\theta^{\prime}\right|_{\infty}$ $\displaystyle=0$ (65)
$\displaystyle\left.\theta^{\prime}\right|_{0}$ $\displaystyle=1$
$\displaystyle\left.\theta^{\prime}\right|_{\infty}$ $\displaystyle=0$ (66)

###  5.1 Scalings

From the previous results we can obtain the corresponding boundary layer
scalings $s$, $t$ and $u$ as those that render equations (59) and (64) and its
boundary conditions (eq. (62) and (65) or (66)) exclusively in terms of
$\eta$, $f$ and $\theta$.

A summary of the scaling factors can be find in the following table:

Author | b | s | t | u
---|---|---|---|---
Blasius[1] | 0 | $1/2$ | $1/2$ | N/A
Ostrach[2] | $\propto C^{-1}\theta$ | $1/4$ | $3/4$ | 0
Sparrow[3] | $\propto C^{-1}\theta x^{u}$ | $1/5$ | $4/5$ | $1/5$
Table 1: <+Caption text+>

###  5.2 Coefficients

Now, we still have the freedom to define the coefficients $A$, $B$ and $C$.
While this is a matter of preference, we report here the most usual practice.
In particular, since we have three unknowns ($A$, $B$ and $C$), we can impose
three new equations, namely for the two terms in which they appear in equation
(LABEL:eqn:self-similiarity-prototype) such that the expression is simplified:

$\displaystyle\frac{1}{\nu AB}$ $\displaystyle=\frac{1}{s}$ (67)
$\displaystyle\frac{bB}{\nu A^{3}}$ $\displaystyle=\theta$ (68)
$\displaystyle\frac{AD}{C}$ $\displaystyle=q$ $\displaystyle C$
$\displaystyle=\frac{1}{c_{wall}}$ (69)

This removes the fractional coefficients in equation (59) and (64)

Author | A | B | C (Dirichlet) | C (Neumann)
---|---|---|---|---
Blasius | $\sqrt{\frac{U}{\nu}}$ | $\frac{1}{\sqrt{\nu U}}$ | $\frac{1}{C_{wall}}$ | $\frac{D}{q}\sqrt{\frac{U}{\nu}}$
Ostrach | $\left(\frac{g\beta\Delta\theta_{w}}{4\nu^{2}}\right)^{1/4}$ | $\left(4^{2}\nu^{2}g\beta\Delta\theta_{w}\frac{q}{D}\right)^{-1/4}$ | $\frac{1}{C_{wall}}$ | N/A
Sparrow | $\left(\frac{1}{5}\frac{g\beta}{\nu^{2}}\frac{q}{D}\right)^{1/5}$ | $\left(5^{4}\nu^{3}g\beta\frac{q}{D}\right)^{-1/5}$ | N/A | $\left(\frac{1}{5}\frac{g\beta}{\nu^{2}}\left(\frac{D}{q}\right)^{4}\right)^{1/5}$
Table 2: Summary of self-similarity scalings and coefficients.

Blasius equation

$\displaystyle\eta=\frac{y}{\delta}=Ayx^{-1/2}$ (70) $\displaystyle f=B\Psi
x^{-1/2}$ (71) $f^{\prime\prime\prime}+2ff^{\prime\prime}=0$ (72)
$\displaystyle f^{\prime}(0)$ $\displaystyle=0$ $\displaystyle f(0)$
$\displaystyle=0$ $\displaystyle f^{\prime}(\infty)$ $\displaystyle=1$ (73)

Ostrach equation

$\displaystyle\eta$ $\displaystyle=\frac{y}{\delta}=Ayx^{-1/4}$ (74)
$\displaystyle f$ $\displaystyle=B\Psi x^{3/4}$ (75) $\displaystyle\theta$
$\displaystyle=C\Delta cx^{-1/4}$ (76) $\displaystyle<+content+>$ (77)

<++>

##  6 Leveque equation

For the advection of passive scalars, when $Sc>>1$ (as it is typically the
case in mass transfer) the mass transfer boundary layer is much thinner than
the hydrodynamic one. This means that a linearization of the velocity profile
is possible as:

$w(x)\approx\frac{\partial w}{\partial x}x=w^{\prime}x$ (78)

This was observed by Leveque, who proposed a new scaling as:

$w^{\prime}x\frac{\partial c}{\partial z}=D\frac{\partial^{2}c}{\partial
x^{2}}$ (79)

where $w^{\prime}$ and $D$ can be functions of $z$ only.

This last assumption is made in order to reuse Blasius-like expressions for
$w^{\prime}$, which then evolves in the streamwise direction $z$.

###  6.1 Dirichlet boundary conditions

To which we can attempt a self-similarity transformation as:

$\eta=Axf(z)^{-s}$ (80)

where now $f(x)$ can be any function of x. On the other hand, the actual
variable $c$ is:

$c=c(\eta)$ (81)

By substituting these into the original equation (79), we obtain:

$\frac{c^{\prime\prime}}{c^{\prime}}=-s\frac{w^{\prime}}{D}f^{\prime}\frac{\eta^{2}}{A^{3}}f^{3s-1}$
(82)

In order to satisfy self-similarity, we impose the following:

$\displaystyle 3s-1$ $\displaystyle=0$ (83)
$\displaystyle\frac{w^{\prime}}{D}f^{\prime}$ $\displaystyle=1$ (84)

from where we obtain the final Leveque scaling:

$\displaystyle s$ $\displaystyle=\frac{1}{3}$ (85) $\displaystyle f$
$\displaystyle=\int\frac{D}{w^{\prime}}dz$ (86)

With this scaling, equation (79) becomes:

$\frac{c^{\prime\prime}}{c^{\prime}}=-\frac{1}{A^{3}}\frac{1}{3}\eta^{2}\\\ $
(87)

which is an eigenvalue problem. We can reduce the order of the problem with
realizing that

$\frac{c^{\prime\prime}}{c^{\prime}}=\frac{\partial}{\partial\eta}ln(c^{\prime})$
(88)

Solving for it we obtain:

$c^{\prime}=Cexp\left(-\frac{1}{A^{3}}\frac{\eta^{3}}{9}\right)$ (89)

From where we can see that, for convenience, we can impose $A$ as:

$A=\left(\frac{1}{9}\right)^{1/3}$ (90)

such that equation (89)

$c^{\prime}=Cexp\left(-\eta^{3}\right)$ (91)

Once $A$ is defined, the similarity is:

$\eta=y\left(\frac{1}{9\int\frac{D}{w^{\prime}}dx}\right)^{\frac{1}{3}}$ (92)

The final expression for $c$ is obtained by integration between $\eta$ and
$\infty$.

$c_{\infty}-c(\eta)=C\int_{\eta}^{\infty}exp\left(-\eta^{3}\right)d\eta=C$
(93)

###  6.2 Neumann boundary conditions

For Neumann boundary conditions:

$\displaystyle D\frac{\partial c}{\partial y}\lvert_{y=0}$
$\displaystyle=\frac{j}{zF}$ (94) $\displaystyle c\lvert_{y\to\infty}$
$\displaystyle=c_{sat}$ (95)

To which we attempt a self-similiarity solution of the form:

$\displaystyle\eta$ $\displaystyle=Ayf(z)^{-s}$ (96)
$\displaystyle\theta(\eta)$ $\displaystyle=B\Delta cf(z)^{-u}$ (97)

where $\Delta c=c-c_{s}at$. We scale $\theta$ as an arbitrary power of the
shape function $f(x)$. Substituting this scaling in the boundary conditions we
obtain:

$\displaystyle f^{u-s}\theta^{\prime}(0)$
$\displaystyle=\frac{j}{zFD}\frac{B}{A}$ (98) $\displaystyle\theta(\infty)$
$\displaystyle=0$ (99)

which requires $u-s=0$ in order to obtain a self-similarity condition. This
means that $\eta$ and $\theta$ scale the same.

Introducing now the scaling in the Leveque equation (79), we obtain the
following ODE:

$w^{\prime}\eta
A^{-1}f^{s+u-1}f^{\prime}\left(u\theta-s\eta\theta^{\prime}\right)=DA^{2}f^{-2s}B^{-1}f^{u}\theta^{\prime\prime}$
(100)

which can be rearranged as

$\frac{w^{\prime}}{D}f^{\prime}f^{3s-1}\eta\left(u\theta-s\eta\theta^{\prime}\right)=A^{3}\theta^{\prime\prime}$
(101)

from where we can impose the following self-similarity conditions:

$\displaystyle\frac{w^{\prime}}{D}f^{\prime}$ $\displaystyle=1$ (102)
$\displaystyle 3s-1$ $\displaystyle=0$ (103)

Self-similarity requires then scaling and the shape function as

$\displaystyle f^{\prime}$ $\displaystyle=\int\frac{D}{w^{\prime}}$ (104)
$\displaystyle s$ $\displaystyle=\frac{1}{3}$ (105)

which fully define the self-similarity. The equation then takes the form

$3A^{3}\theta^{\prime\prime}+\eta^{2}\theta^{\prime}-\eta\theta=0$ (106)

where $A$ can be chosen to simplify the system as

$A=\left(\frac{1}{9}\right)^{1/3}$ (107)

which renders the final form of the Leveque equation with Neumann boundary
conditions as:

$\frac{1}{3}\theta^{\prime\prime}+\eta^{2}\theta^{\prime}-\eta\theta=0$ (108)

which has an analytical solution if we attempt a solution like:

$\theta(\eta)=\eta v(\eta)$ (109)

which accepts an analytical solution of the form

$\theta(\lambda)=CE_{\frac{4}{3}}\left(\lambda^{3}\right)=C\lambda\Gamma\left(-\frac{1}{3},\lambda^{3}\right)$
(110)

where $E_{i}$ is the exponential integral, while $\Gamma$ is the incomplete
gamma function. Taking its boundary conditions we obtain:

$\displaystyle\theta^{\prime}(0)$
$\displaystyle=\Gamma\left(-\frac{1}{3}\right)C$ (111)
$\displaystyle\theta(\infty)$ $\displaystyle=0$ (112)

Thus, the final form of $\theta$ becomes:

$\theta(\eta)=\eta\frac{\Gamma\left(-\frac{1}{3},x^{3}\right)}{\Gamma\left(-\frac{1}{3}\right)}=\eta
Q\left(-\frac{1}{3},x^{3}\right)$ (113)

where $Q$ is the regularized upper incomplete gamma function.

The final expression for mass then becomes:

$\Delta
c_{wall}=\frac{q}{D}\theta(\eta)\left(9\int_{0}^{z}\frac{D}{w^{\prime}}dt\right)^{1/3}$
(114)

which yields the expression for mass transfer coefficient as:

$k_{f}=\frac{1}{\theta(0)}\left(9\int_{0}^{z}\frac{D}{w^{\prime}}dt\right)^{-1/3}$
(115)

Where the final Sherwood number is:

$Sh_{f}=\frac{k_{f}z}{D}=\frac{1}{\theta(0)}z\left(9\int_{0}^{z}\frac{D}{w^{\prime}}dt\right)^{-1/3}$
(116)

## References

  * [1] H. Blasius.  The boundary layer in fluids with little friction (J. Vanier, Trans.).  Technical report, National Advisory Committee for Aeronautics, Washington, 1950.
  * [2] Simon Ostrach.  An analysis of laminar free-convection flow and heat transfer about a flat plate parallel to the direction of the generating body force.  Technical report, National Advisory Commitee for Aeronautics, 1953.
  * [3] E. M. Sparrow, R. Eichhorn, and J. L. Gregg.  Combined forced and free convection in a boundary layer flow.  Phys. Fluids, 2(3):319–328, 1959.