# Monte-Carlo method for incompressible fluid
flows past obstacles

By V. Cherepanov and Z. Qian Mathematical Institute, University of Oxford,
Oxford OX2 6GG. Email<EMAIL_ADDRESS>Institute, University of Oxford, Oxford OX2 6GG and Oxford Suzhou Centre for
Advanced Research, Suzhou, China. Email<EMAIL_ADDRESS>
###### Abstract

We establish stochastic functional integral representations for incompressible
fluid flows occupying wall-bounded domains using the conditional law duality
for a class of diffusion processes. These representations are used to derive a
Monte-Carlo scheme based on the corresponding exact random vortex formulation.
We implement several numerical experiments based on the Monte-Carlo method
without appealing to the boundary layer flow computations, to demonstrate the
methodology.

_Key words_ : diffusion processes, incompressible fluid flow, Monte-Carlo
simulation, random vortex method

_MSC classifications_ : 76M35, 76M23, 60H30, 65C05, 68Q10.

##  1 Introduction

In this paper we aim to develop Monte-Carlo schemes for the numerical analysis
of incompressible flows past solid obstacles, that is flows occupying certain
domains with a nontrivial boundary. These fluid flows satisfy the
incompressible Navier-Stokes equations, and we derive exact functional
integral representations for their solutions. Phrasing these in terms of
random vortex formulations, we are able to establish corresponding numerical
schemes and conduct simulations for such flows.

The flows in question are characterised by their velocity $u(x,t)$ for $x\in
D$, where $D$ is a domain in $\mathbb{R}^{2}$ or $\mathbb{R}^{3}$ with
nontrivial boundary $\partial D$. The velocity then follows the Navier-Stokes
equations:

$\frac{\partial}{\partial t}u+(u\cdot\nabla)u-\nu\Delta u+\nabla P=F,$ (1.1)
$\nabla\cdot u=0$ (1.2)

in $D$, and the no-slip condition for $u$ is satisfied, i.e. $u(x,t)=0$ for
$x\in\partial D$. The constant $\nu>0$ is the kinematic viscosity of the
fluid, $P$ is the pressure and $F$ denotes an external force.

One can get another important formulation of the problem, writing it in terms
of different dynamical variables, namely the velocity $u$ and the vorticity
$\omega=\nabla\wedge u$, the curl of the velocity. For simplicity, we assume
in the following derivation that the domain $D$ is two-dimensional — we refer
the reader to Section 2 for the three-dimensional case. We also assume that
the velocity $u$ is extended to the whole space $\mathbb{R}^{2}$ such that
$\nabla\cdot u=0$ in distribution (e.g., by letting $u=0$ outside of the
domain $D$, see Sections 2 and 4 for details).

Indeed, the equation (1.1) implies the vorticity transport equation

$\frac{\partial}{\partial t}\omega+(u\cdot\nabla)\omega-\nu\Delta\omega=G$
(1.3)

in $D$, where $G=\nabla\wedge F$. Notice, however, that in general the value
of the vorticity $\omega$ is non-zero along the boundary $\partial D$, and we
denote $\theta=\left.\omega\right|_{\partial D}$. It is well-known that the
solution to equations of the form (1.3) can be expressed in terms of
functional integrals — this idea may trace back to Feynman [13] and Kac [25].
However, we first have to introduce a perturbation $W_{\varepsilon}$ of the
vorticity $\omega$ such that $W_{\varepsilon}=\omega$ everywhere except a thin
$\varepsilon$-layer of the boundary $\partial D$, and
$\left.W_{\varepsilon}\right|_{\partial D}=0$. One way of defining such
perturbation is given in [38] using a smooth cutoff function, which we also
utilise in the subsequent sections.

Therefore, for $W_{\varepsilon}$ the vorticity equation (1.3) implies

$\frac{\partial}{\partial
t}W_{\varepsilon}+(u\cdot\nabla)W_{\varepsilon}-\nu\Delta
W_{\varepsilon}=g_{\varepsilon}$ (1.4)

in $D$ where $\left.W_{\varepsilon}\right|_{\partial D}=0$ and
$g_{\varepsilon}$ is the corresponding perturbation of $G$. Introducing the
Green function $p(s,x,t,y)$ for the forward parabolic operator
$\nu\Delta-u\cdot\nabla-\frac{\partial}{\partial t}$ in $D$ subject to the
Dirichlet boundary condition, we write

$W_{\varepsilon}(y,t)=\int_{D}p(0,\xi,t,y)W_{\varepsilon}(\xi,0)\textrm{d}\xi+\int_{0}^{t}\int_{D}p(s,\xi,t,y)g_{\varepsilon}(\xi,s)\textrm{d}\xi\textrm{d}s$
(1.5)

which we can then use to recover the unperturbed vorticity $\omega$. Note that
for incompressible flows one then recovers the velocity by

$u(x,t)=\int_{D}K(x,y)\omega(y,t)\textrm{d}y$ (1.6)

where $K$ is the Biot-Savart kernel of the domain $D$ (assuming both the
velocity and the vorticity decay to zero at infinity sufficiently fast).

To give a formulatic expression, we assume throughout the introduction that a
representation similar to (1.5) holds for the unperturbed vorticity $\omega$,
that is

$\omega(y,t)=\int_{D}p(0,\xi,t,y)\omega(\xi,0)\textrm{d}\xi+\int_{0}^{t}\int_{D}p(s,\xi,t,y)g(\xi,s)\textrm{d}\xi\textrm{d}s.$
(1.7)

Then for the velocity, one has

$u(x,t)=\int_{D}\int_{D}K(x,y)p(0,\xi,t,y)\omega(\xi,0)\textrm{d}\xi\textrm{d}y+\int_{0}^{t}\int_{D}\int_{D}K(x,y)p(s,\xi,t,y)g(\xi,s)\textrm{d}\xi\textrm{d}s\textrm{d}y.$
(1.8)

We stress again that even though the above formulae hold for the case of
unbounded flows, in the case of nontrivial boundary the actual representations
essentially involve the boundary vorticity $\theta$, see Sections 2 and 4 for
exact expressions.

It turns out one can formulate the above as stochastic representations due to
the following observation. Noting again that the vector field $u$ is
divergence-free, we have that the function $p(s,x,t,y)$ is also the transition
density for the diffusion process with infinitesimal generator
$\nu\Delta+u\cdot\nabla$. This means that for the family of stochastic
processes $X^{\xi,s}$ given by

$\textrm{d}X_{t}^{\xi,s}=u(X^{\xi,s}_{t},t)\textrm{d}t+\sqrt{2\nu}\textrm{d}B_{t}$
(1.9)

for $t\geq s$ and $X^{\xi,s}_{t}=\xi$ for $t\leq s$, we have that

$p(s,\xi,t,y)\textrm{d}y=\mathbb{P}\left[X^{\xi,s}_{t}\in\textrm{d}y\right].$
(1.10)

The processes are called Taylor’s diffusions (due to Taylor [43]), one
interprets the processes as "imaginary" Brownian particles following the flow
with velocity $u(x,t)$.

Let us view the integral representations we have so far in terms of
expectations with respect to the family of Taylor’s diffusions introduced
above. Indeed, it is easy to see that the integrals in the representation
(1.8) can be viewed in terms of the expectations with respect to $X^{\xi,s}$,
that is

$u(x,t)=\int_{D}\mathbb{E}\left[K(x,X^{\xi,0}_{t})\right]\omega(\xi,0)\textrm{d}\xi+\int_{0}^{t}\int_{D}\mathbb{E}\left[K(x,X^{\xi,s}_{t})\right]g(\xi,s)\textrm{d}\xi\textrm{d}s.$
(1.11)

Having discussed the above representations, we propose the following numerical
scheme. We couple the processes $X^{\xi,s}$ with the velocity $u(x,t)$ and
update them according to approximations given by (1.9) and (1.11)
respectively. Indeed, the former depends only on the velocity $u(x,t)$ for
which one can take the approximation given by the latter. Note that, in turn,
the velocity is determined by the diffusions $X^{\xi,s}$ and the presence of
expectations in the representation (1.11) allows for Monte-Carlo-type schemes
(we refer the reader to the Section 5 for particular schemes we used in
simulations).

One of the drawbacks for the scheme described above is that it involves
Taylor’s diffusions (1.9) started at every time $s\in[0,t]$. This implies that
the memory required for the computation and the computation complexity itself
are increasing with time $t$. In Section 2, we overcome this difficulty by
deriving representations different from (1.7) and (1.8) using the duality of
the conditional laws of Taylor’s diffusions proved in [37]. Although the
representations in this case involve more complicated integrands (see Theorems
2.1, 2.2 and their two-dimensional versions, Theorems 2.3, 2.4), we are able
to devise the schemes with non-increasing memory and computational complexity
in Section 5.

We also have to mention again that the above representations (1.7) and (1.8)
in the stated form hold for domains without boundary. In the case of domains
with nontrivial boundary, the representations are derived from (1.5) for the
perturbed vorticity $W_{\varepsilon}$. These depend on the boundary vorticity
$\theta$ which poses an additional difficulty in simulations as we are unable
to compute the values of $\theta$ using the expressions for $\omega$. However,
as in practice we use a desingularised kernel $K_{\delta}$ to compute the
velocity, we write the formal derivative of the representation for $u(x,t)$ to
compute the vorticity $\omega$ used in plots and the boundary vorticity
$\theta$.

Note also that as the expression (1.5) for the perturbed vorticity
$W_{\varepsilon}$ depends on the boundary layer thickness parameter
$\varepsilon$, it is interesting to consider the limit of the representations
derived from (1.5) as $\varepsilon\to 0$. We do the corresponding computations
in Section 4 for different domains and notice that many terms in fact do not
contribute to the limit. This observation turns out to be useful in Section 5
as it simplifies the numerical schemes significantly.

Let us remark that the discussed above limiting representations obtained in
Section 4 have a certain interest on their own as they display how the
boundary vorticity $\theta$ influences the velocity $u$ (see Section 4 for
details), which seems to be in accordance with Prandtl’s boundary layer theory
[35]. In some of the experiments reported in Section 5 we omit the terms
involving $\theta$ when we compute the velocity $u$ and notice that they
essentially contribute to the chaotic behaviour of the flow close to the
boundary.

We also note that we derive these representations in Section 4 for two
particular domains (namely, the ones for flows passing a flat plate and a
wedge) and for general two-dimensional domains $D$ conformally equivalent to
the half-plane with smooth boundary $\partial D$. However, we believe that
similar representations should hold for more complicated domains as well
(e.g., not necessary simply connected domains with piecewise smooth boundary).

For the existing literature on the subject, we first mention that there has
been extensive literature studying Navier-Stokes equations, in particular,
numerical analysis of their solutions. For the general numerical methods, we
mention Computational Fluid Dynamics (CFD) (see [19], [45]), while more
particular numerical approaches include Direct Numerical Simulations (DNS)
(see [33], [41], [32] and also [45]), Large Eddy Simulations (LES) (see [15],
[29], [27]), Probability Density Function (PDF) (see [34]). There is also a
large volume of literature with numerical study of turbulence, e.g., see the
works [3], [22], [14], [23], [39], [42], [47, 48, 49], as well as those that
study other various aspects of boundary flows [2], [16], [17], [18], [24],
[26], [40], [46].

Random vortex method is based originally on the work [4] by Chorin, though the
idea of using Brownian fluid particles dates back to Taylor [43]. It is also
successfully used in studying different aspects of turbulent flows, see [34],
[31] and [12]. We also mention that some convergence results are known for
random vortex methods, see [1], [21] and [30], however, these concern the
flows without boundary and the convergence results for the schemes described
above will be studied in a future work. For other probabilistic aspects of
fluid dynamics we refer the reader to [5, 6] and [7], and we also mention
particularly stochastic Lagrangian approach that is used to study isotropic
turbulence, see [8, 9], [10, 11].

This paper is based on the previous works [36], [37] and [38] and generalises
the results for wall-bounded domains presented in [36] and [38]. In the recent
paper [28], the ideas similar to those presented here are used to study
Oberbeck-Boussinesq flows which are outside of the scope of the current work.

The paper is organised as follows. In Section 2, we derive a random vortex
representation for general two- and three-dimensional flows occupying the
half-space. This representation is derived using the duality of conditional
laws of the Taylor diffusions and thus depends only on diffusions started at
time $0$. In Section 3, we review the Biot-Savart law and state auxiliary
results concerning the Biot-Savart kernel for the domains in question. In
Section 4, we consider representations as those described above and compute
their limits as the thickness of the thin layers converges to zero. In Section
5, we use the above results to derive numerical schemes and provide the
experiment results.

##  2 Random vortex for wall-bounded flows

In this section, we follow the approach proposed in [36] and [38]. In the
latter paper the random vortex dynamics has been established for viscous fluid
flows moving along with a solid wall modelled by the two dimensional Navier-
Stokes equations. Let $u=(u^{1},u^{2},u^{3})$ denote the velocity of an
incompressible viscous fluid flow constrained in the upper half space $D$
where $x_{3}\geq 0$. The velocity $u$ satisfies the no-slip condition, i.e.
$u(x,t)=0$ for $x\in\partial D$ and $t>0$. Let us assume that there is an
external force $F=(F^{1},F^{2},F^{3})$ supplying the energy to the fluid
dynamical system. Therefore the velocity $u$ and the pressure $P$ are evolved
according to the Navier-Stokes equations:

$\frac{\partial u^{i}}{\partial t}+\sum_{j=1}^{3}u^{j}\frac{\partial
u^{i}}{\partial x_{j}}-\nu\Delta u^{i}+\frac{\partial P}{\partial
x_{i}}-F^{i}=0\quad\textrm{ in }D$ (2.1)

for $i=1,2,3$, and

$\sum_{j=1}^{3}\frac{\partial u^{j}}{\partial x_{j}}=0\quad\textrm{ in }D.$
(2.2)

The initial velocity is denoted by $u_{0}(x)=u(x,0)$. In the following we will
often omit the summation over repeated indices.

To proceed our discussion, we assume that the velocity $u$ is smooth inside
$D$ and is $C^{3}$ up to the boundary $\partial D$. This assumption is
although technical, but in no means it is trivial and apparent. Indeed it
remains and will still be stand as an open problem in mathematics for
establishing the regularity for solutions to the general three dimensional
Navier-Stokes equations. We do not pursue this line of research which lies
outside the scope of current project.

We introduce the vorticity $\omega=\nabla\wedge u$, whose components
$\omega^{i}=\varepsilon^{ijk}\frac{\partial}{\partial x^{j}}u^{k}$. In our
study the vorticity transport equations play a crucial role. The vorticity
transport equations are evolution equations for $\omega$, which are given as
the following

$\frac{\partial\omega^{i}}{\partial t}+u^{j}\frac{\partial\omega^{i}}{\partial
x_{j}}-\nu\Delta\omega^{i}+\omega^{j}S_{j}^{i}-G^{i}=0\quad\textrm{ in }D$
(2.3)

where $S_{j}^{i}=\frac{1}{2}(\frac{\partial u^{i}}{\partial
x_{j}}+\frac{\partial u^{j}}{\partial x_{i}})$ is the symmetric tensor of
rate-of-strain, and $G=\nabla\wedge F$.

We adopt the simple idea in [38] of extending the definition of the velocity
$u=(u^{1},u^{2},u^{3})$ to the whole space $\mathbb{R}^{3}$ by the reflection
principle, that is

$u^{i}(\bar{x},t)=u^{i}(x,t)\quad\textrm{ for }i=1,2;\quad
u^{3}(\bar{x},t)=-u^{3}(x,t)$

for $x_{3}<0$. Since $u$ satisfies the no-slip condition, this yields a
crucial fact that the extended velocity $u$ is divergence free in the
distribution sense on $\mathbb{R}^{3}$. Therefore the adjoint of the heat
operator $\nu\Delta-u\cdot\nabla-\frac{\partial}{\partial t}$ is just the heat
operator $\nu\Delta+u\cdot\nabla+\frac{\partial}{\partial t}$. Therefore the
Green function $\Gamma_{D}(\tau,x,t,y)$ to the Dirichlet boundary problem of
the (forward) parabolic equation

$\begin{cases}\left(\frac{\partial}{\partial
t}-\nu\Delta+u\cdot\nabla\right)f=0&\textrm{ in }D,\\\ f(x,t)=0&\textrm{ if
}x_{3}=0\end{cases}$

coincides with the transition probability density of the diffusion process
with its infinitesimal generator $\nu\Delta+u\cdot\nabla$ killed on leaving
the domain $D$. This leads to the following construction. Let $X$ be the
diffusion process with infinitesimal generator $\nu\Delta+u\cdot\nabla$, which
is a diffusion with state space $\mathbb{R}^{3}$, without “killing” at the
boundary $x_{3}=0$. That is, $X$ is a weak solution of the stochastic
differential equation

$\textrm{d}X=u(X,t)\textrm{d}t+\sqrt{2\nu}\textrm{d}B.$ (2.4)

Let $p(\tau,x,t,y)$ (for $t>\tau\geq 0$, and $x,y\in\mathbb{R}^{3}$) be the
transition probability density function, which is positive and Hölder’s
continuous in all arguments, as long as $u$ is bounded and Borel measurable.
Formally

$p(\tau,x,t,y)\textrm{d}y=\mathbb{P}\left[X_{t}\in\textrm{d}y|X_{\tau}=x\right]$

and therefore it is clear that

$p(\tau,x,t,y)=p(\tau,\overline{x},t,\overline{y})\quad\textrm{ for
}x,y\in\mathbb{R}^{3}.$

Hence, by applying the reflection principle we have therefore the following
representation

$\Gamma_{D}(\tau,x,t,y)=p(\tau,x,t,y)-p(\tau,x,t,\overline{y})\quad\textrm{
for }x,y\in D.$ (2.5)

We next apply this representation to the study of the vorticity $\omega$. The
vorticity $\omega$ may be considered as a solution of the linear parabolic
equation (2.3) if $u$ is supposed as a given fluid dynamic variable. However,
$\omega$ has nontrivial boundary value in general, so let $\theta$ be the
trace of $\omega$ along the boundary $\partial D$ (note that the vorticity
boundary values may be identified with the stress of the fluid flow
immediately injected to the wall, cf. [40] [44]). Then it is easy to see that
$\omega^{3}$ has trace zero, and due to the no-slip condition imposed on $u$,
$\theta$ can be identified with the normal part of the stress applied to the
boundary $\partial D$. Note that $\theta$ is a time dependent vector field on
the boundary $\partial D$.

To handle the no vanishing boundary vorticity $\theta$, we employ the same
technique used in [38]. Let $\phi$ be a smooth cut-off function defined on
$[0,\infty)$ with values in $[0,1]$, such that $\phi(r)=1$ for $r\in[0,1/3]$
and $\phi(r)=0$ for $r\geq 2/3$, and define
$\sigma_{\varepsilon}(x,t)=\phi(x_{3}/\varepsilon)\theta(x_{1},x_{2},t)$ for
$x=(x_{1},x_{2},x_{3})\in D$. Let

$W_{\varepsilon}(x,t)=\omega(x,t)-\sigma_{\varepsilon}(x,t).$ (2.6)

Then

$\frac{\partial W_{\varepsilon}^{i}}{\partial t}+u^{j}\frac{\partial
W_{\varepsilon}^{i}}{\partial x_{j}}-\nu\Delta
W_{\varepsilon}^{i}-W_{\varepsilon}^{j}\frac{\partial u^{i}}{\partial
x_{j}}-g_{\varepsilon}^{i}=0\quad\textrm{ in }D,$ (2.7)

for $i=1,2,3$, and

$\left.W_{\varepsilon}(x,t)\right|_{x\in\partial D}=0\textrm{ }\quad\textrm{
for }x\in\partial D,$ (2.8)

where

$g_{\varepsilon}^{i}=G^{i}-\frac{\partial}{\partial
t}\sigma_{\varepsilon}^{i}-u^{j}\frac{\partial}{\partial
x_{j}}\sigma_{\varepsilon}^{i}+\nu\Delta\sigma_{\varepsilon}^{i}-\sigma_{\varepsilon}^{j}S_{j}^{i}$
(2.9)

for $i=1,2,3$.

###  2.1 Random vortex for wall-bounded flows

According to the Feynman-Kac formula for forward heat equations established in
[36] we can represent $W_{\varepsilon}$ in terms of the distribution of the
Taylor diffusion. Here we present a slightly different approach. To this end,
we have to introduce some notation. Let $\mathbb{P}^{\xi}$ denote the
distribution of the diffusion process with infinitesimal generator
$\nu\Delta+u\cdot\nabla$ (where $u$ is extended on $\mathbb{R}^{3}$ by the
reflection principle) started from $\xi\in\mathbb{R}^{3}$ at time $0$, and
$\mathbb{P}_{t}^{\xi\rightarrow\eta}$ the conditional law of
$\mathbb{P}^{\xi}\left[\cdot|\psi(t)=\eta\right]$ where $\psi$ is the
coordinate process (i.e. canonical element) on the path space
$C([0,\infty),\mathbb{R}^{3})$.

We recall that $p_{b}(s,x,t,y)$ denotes the transition probability density
function of the diffusion with generator $\nu\Delta+b\cdot\nabla$.

######  Theorem 2.1.

Let $X^{\eta}$ be the Taylor diffusion:

$\textrm{d}X_{t}^{\eta}=u(X_{t}^{\eta},t)\textrm{d}t+\sqrt{2\nu}\textrm{d}B_{t},\quad
X_{0}^{\eta}=\eta$ (2.10)

for every $\eta\in\mathbb{R}^{3}$. For each pair $\eta\in\mathbb{R}^{3}$ and
$t>0$, define $s\mapsto Q(\eta,t;s)$ to be the unique solution of the
differential equations:

$\frac{\textrm{d}}{\textrm{d}s}Q_{j}^{i}(\eta,t;s)=-Q_{k}^{i}(\eta,t;s)1_{D}(X_{s}^{\eta})A_{j}^{k}(X_{s}^{\eta},s),\quad
Q_{j}^{i}(\eta,t;t)=\delta_{ij}$ (2.11)

where $A_{j}^{k}=\frac{\partial}{\partial x_{j}}u^{k}$ and $i,j,k=1,2,3$. Then

$\displaystyle\omega^{i}(\xi,t)$
$\displaystyle=\sigma_{\varepsilon}^{i}(\xi,t)+\int_{D}\mathbb{P}^{\eta\rightarrow\xi}\left[Q_{j}^{i}(\eta,t;0)1_{\\{t<\zeta(X^{\eta}\circ\tau_{t})\\}}\right]W_{\varepsilon}^{j}(\eta,0)p_{u}(0,\eta,t,\xi)\textrm{d}\eta$
$\displaystyle+\int_{0}^{t}\int_{D}\mathbb{P}^{\eta\rightarrow\xi}\left[Q_{j}^{i}(\eta,t;s)1_{\\{t-s<\zeta(X^{\eta}\circ\tau_{t})\\}}g_{\varepsilon}^{j}(X_{s}^{\eta},s)\right]p_{u}(0,\eta,t,\xi)\textrm{d}\eta\textrm{d}s$
(2.12)

for every $\xi\in D$ and $t>0$, $i=1,2,3$, where
$\zeta(\psi)=\inf\left\\{s:\psi(s)\notin D\right\\}$ and $\tau_{t}$ denotes
the time reversion operator on the path space $C([0,t];\mathbb{R}^{3})$.

###### Proof.

For simplicity, let $q_{j}^{i}(x,t)=1_{D}(x)A_{j}^{i}(x,t)$, where
$i,j=1,2,3$. Let $T>0$ be fixed. Recall that $u(x,t)$ is extended for all
$x\in\mathbb{R}^{3}$ such that $\nabla\cdot u=0$ in distribution on
$\mathbb{R}^{3}$. Let $\tilde{X}^{\xi}$ be the solution to the stochastic
differential equation

$\textrm{d}\tilde{X}_{t}^{\xi}=-u(\tilde{X}_{t}^{\xi},T-t)\textrm{d}t+\sqrt{2\nu}\textrm{d}B_{t},\quad\tilde{X}_{0}^{\xi}=\xi$
(2.13)

which is understood as the stochastic integral equation

$\tilde{X}_{t}^{\xi}=\xi-\int_{0}^{t\wedge
T}u(\tilde{X}_{s}^{\xi},T-s)\textrm{d}s+\sqrt{2\nu}\int_{0}^{t\wedge
T}\textrm{d}B_{s}\quad\textrm{ for all }t\geq 0.$ (2.14)

Let

$T_{\xi}=\inf\left\\{t\geq 0:\tilde{X}_{t}^{\xi}\notin D\right\\}$

be the first time the process $\tilde{X}^{\xi}$ leaves the region $D$. Then
$t\rightarrow\tilde{X}_{t\wedge T_{\xi}}^{\xi}$ is a diffusion process too,
and

$\displaystyle\tilde{X}_{t\wedge T_{\xi}}^{\xi}$
$\displaystyle=\xi-\int_{0}^{t\wedge T\wedge
T_{\xi}}u(\tilde{X}_{s}^{\xi},T-s)\textrm{d}s+\sqrt{2\nu}\int_{0}^{t\wedge
T\wedge T_{\xi}}\textrm{d}B_{s}$
$\displaystyle=\xi-\int_{0}^{t}1_{\left\\{s<T\wedge
T_{\xi}\right\\}}u(\tilde{X}_{s}^{\xi},T-s)\textrm{d}s+\sqrt{2\nu}\int_{0}^{t}1_{\left\\{s<T\wedge
T_{\xi}\right\\}}\textrm{d}B_{s}$ (2.15)

for all $t\geq 0$. According to Feynman-Kac, we define

$\textrm{d}\tilde{Q}_{j}^{i}(s)=\tilde{Q}_{k}^{i}(s)q_{j}^{k}(\tilde{X}_{s}^{\xi},T-s)\textrm{d}s,\quad\tilde{Q}_{j}^{i}(0)=\delta_{ij}$
(2.16)

where $i,j=1,2,3$. Let $Y_{t}^{j}=W_{\varepsilon}^{j}(\tilde{X}_{t\wedge
T_{\xi}}^{\xi},T-t)$ which is well defined, as according to our assumption,
$W_{\varepsilon}$ is $C^{2}$ on $\overline{D}$. Since $W_{\varepsilon}$
vanishes along the boundary $\partial D$, so that

$\displaystyle Y_{t}^{j}$
$\displaystyle=W_{\varepsilon}^{j}(\tilde{X}_{t\wedge
T_{\xi}}^{\xi},T-t)=1_{\\{t<T_{\xi}\\}}W_{\varepsilon}^{j}(\tilde{X}_{t}^{\xi},T-t)+1_{\\{t\geq
T_{\xi}\\}}W_{\varepsilon}^{j}(\tilde{X}_{T_{\xi}}^{\xi},T-t)$
$\displaystyle=1_{\\{t<T_{\xi}\\}}W_{\varepsilon}^{j}(\tilde{X}_{t}^{\xi},T-t)$

for all $t\leq T$. Moreover, according to Itô’s formula,

$\textrm{d}Y_{t}^{j}=\nabla W_{\varepsilon}^{j}(\tilde{X}_{t\wedge
T_{\xi}}^{\xi},T-t)\cdot\textrm{d}\tilde{X}_{t\wedge
T_{\xi}}^{\xi}-\frac{\partial W_{\varepsilon}^{j}}{\partial
t}(\tilde{X}_{t\wedge T_{\xi}}^{\xi},T-t)\textrm{d}t+\nu\Delta
W_{\varepsilon}^{j}(\tilde{X}_{t\wedge T_{\xi}}^{\xi},T-t)\textrm{d}t.$

Since $W_{\varepsilon}$ vanishes identically on $\partial D$,

$\displaystyle\frac{\partial W_{\varepsilon}^{j}}{\partial
t}(\tilde{X}_{t\wedge T_{\xi}}^{\xi},T-t)$
$\displaystyle=1_{\left\\{t<T_{\xi}\right\\}}\frac{\partial
W_{\varepsilon}^{j}}{\partial t}(\tilde{X}_{t}^{\xi},T-t)+1_{\left\\{t\geq
T_{\xi}\right\\}}\frac{\partial W_{\varepsilon}^{j}}{\partial
t}(\tilde{X}_{T_{\xi}}^{\xi},T-t)$
$\displaystyle=1_{\left\\{t<T_{\xi}\right\\}}\frac{\partial
W_{\varepsilon}^{j}}{\partial t}(\tilde{X}_{t}^{\xi},T-t)$

and similarly

$\displaystyle\nabla W_{\varepsilon}^{j}(\tilde{X}_{t\wedge
T_{\xi}}^{\xi},T-t)$ $\displaystyle=1_{\left\\{t<T_{\xi}\right\\}}\nabla
W_{\varepsilon}^{j}(\tilde{X}_{t}^{\xi},T-t),$ $\displaystyle\Delta
W_{\varepsilon}^{j}(\tilde{X}_{t\wedge T_{\xi}}^{\xi},T-t)$
$\displaystyle=1_{\left\\{t<T_{\xi}\right\\}}\Delta
W_{\varepsilon}^{j}(\tilde{X}_{t}^{\xi},T-t)$

for all $t\leq T$. Therefore

$\displaystyle\textrm{d}Y_{t}^{j}$
$\displaystyle=\sqrt{2\nu}1_{\left\\{t<T_{\xi}\right\\}}\nabla
W_{\varepsilon}^{j}(\tilde{X}_{t}^{\xi},T-t)\cdot\textrm{d}B_{t}$
$\displaystyle+1_{\left\\{t<T_{\xi}\right\\}}\left(\nu\Delta
W_{\varepsilon}^{j}-(u\cdot\nabla W_{\varepsilon}^{j})-\frac{\partial
W_{\varepsilon}^{j}}{\partial t}\right)(\tilde{X}_{t}^{\xi},T-t)\textrm{d}t$
(2.17)

for all $t\leq T$. By (2.7) which may be written as the following

$\left(\nu\Delta-u\cdot\nabla-\frac{\partial}{\partial
t}\right)W_{\varepsilon}^{j}=-A_{k}^{j}W_{\varepsilon}^{k}-g_{\varepsilon}^{j}\quad\textrm{
in }D$ (2.18)

hence, by substituting this into (2.17), it follows that

$\textrm{d}Y_{t}^{j}=\sqrt{2\nu}1_{\\{t<T_{\xi}\\}}\nabla
W_{\varepsilon}^{j}(\tilde{X}_{t}^{\xi},T-t)\cdot\textrm{d}B_{t}-1_{\\{t<T_{\xi}\\}}\left(A_{k}^{j}W_{\varepsilon}^{k}+g_{\varepsilon}^{j}\right)(\tilde{X}_{t}^{\xi},T-t)\textrm{d}t.$
(2.19)

Let $M^{i}=\tilde{Q}_{j}^{i}Y^{j}$. Using stochastic integration by parts one
obtains

$\displaystyle M_{t}^{i}$
$\displaystyle=\tilde{Q}_{j}^{i}(0)Y_{0}^{j}+\int_{0}^{t}\tilde{Q}_{j}^{i}(s)\textrm{d}Y_{s}^{j}+\int_{0}^{t}Y_{s}^{j}\textrm{d}\tilde{Q}_{j}^{i}(s)$
$\displaystyle=\tilde{Q}_{j}^{i}(0)Y_{0}^{j}+\sqrt{2\nu}\int_{0}^{t}\tilde{Q}_{j}^{i}(s)1_{\\{s<T_{\xi}\\}}\nabla
W_{\varepsilon}^{j}(\tilde{X}_{s}^{\xi},T-s)\cdot\textrm{d}B_{s}$
$\displaystyle-\int_{0}^{t}\tilde{Q}_{j}^{i}(s)1_{\\{s<T_{\xi}\\}}\left(A_{k}^{j}W_{\varepsilon}^{k}+g_{\varepsilon}^{j}\right)(\tilde{X}_{s}^{\xi},T-s)\textrm{d}s+\int_{0}^{t}W_{\varepsilon}^{j}(\tilde{X}_{s\wedge
T_{\xi}}^{\xi},T-s)\textrm{d}\tilde{Q}_{j}^{i}(s)$
$\displaystyle=\tilde{Q}_{j}^{i}(0)Y_{0}^{j}+\sqrt{2\nu}\int_{0}^{t}\tilde{Q}_{j}^{i}(s)1_{\\{s<T_{\xi}\\}}\nabla
W_{\varepsilon}^{j}(\tilde{X}_{s}^{\xi},T-s)\cdot\textrm{d}B_{s}$
$\displaystyle-\int_{0}^{t}1_{\\{s<T_{\xi}\\}}\tilde{Q}_{j}^{i}(s)g_{\varepsilon}^{j}(\tilde{X}_{s}^{\xi},T-s)\textrm{d}s$
$\displaystyle+\int_{0}^{t}1_{\\{s<T_{\xi}\\}}W_{\varepsilon}^{j}(\tilde{X}_{s}^{\xi},T-s)\left(-\tilde{Q}_{k}^{i}(s)q_{j}^{k}(\tilde{X}_{s}^{\xi},T-s)\textrm{d}s+\textrm{d}\tilde{Q}_{j}^{i}(s)\right)$

here the last equality follows from the fact that on
$\left\\{s<T_{\xi}\right\\}$, the process $\tilde{X}_{s}^{\xi}$ takes values
in $D$, so that $A_{j}^{i}$ coincides with $q_{j}^{i}$. This computation holds
for any differentiable, adapted processes $\tilde{Q}_{j}^{i}$. In particular,
if $\tilde{Q}$ solves (2.16), then

$\displaystyle M_{t}^{i}$
$\displaystyle=Y_{0}^{i}+\sqrt{2\nu}\int_{0}^{t}1_{\\{s<T_{\xi}\\}}\tilde{Q}_{j}^{i}(s)\nabla
W_{\varepsilon}^{j}(\tilde{X}_{s}^{\xi},T-s)\cdot\textrm{d}B_{s}$
$\displaystyle-\int_{0}^{t}1_{\\{s<T_{\xi}\\}}\tilde{Q}_{j}^{i}(s)g_{\varepsilon}^{j}(\tilde{X}_{s}^{\xi},T-s)\textrm{d}s$

which yields that

$\displaystyle M_{T\wedge T_{\xi}}^{i}$
$\displaystyle=W_{\varepsilon}^{i}(\xi,T)-\int_{0}^{T\wedge
T_{\xi}}1_{\\{t<T_{\xi}\\}}\tilde{Q}_{j}^{i}(t)g_{\varepsilon}^{j}(\tilde{X}_{t}^{\xi},T-t)\textrm{d}t$
$\displaystyle+\sqrt{2\nu}\int_{0}^{T\wedge
T_{\xi}}1_{\\{t<T_{\xi}\\}}\tilde{Q}_{j}^{i}(t)\nabla
W_{\varepsilon}^{j}(\tilde{X}_{t}^{\xi},T-t)\cdot\textrm{d}B.$ (2.20)

Since

$M_{T\wedge
T_{\xi}}^{i}=1_{\\{T<T_{\xi}\\}}\tilde{Q}_{j}^{i}(T)W_{\varepsilon}^{j}(\tilde{X}_{T}^{\xi},0),\quad
Y_{0}^{i}=W_{\varepsilon}^{i}(\xi,T)$

so that, after taking expectation of the both sides of (2.20) to obtain that

$W_{\varepsilon}^{i}(\xi,T)=\mathbb{E}\left[1_{\\{T<T_{\xi}\\}}\tilde{Q}_{j}^{i}(T)W_{\varepsilon}^{j}(\tilde{X}_{T}^{\xi},0)\right]+\mathbb{E}\left[\int_{0}^{T}1_{\\{t<T_{\xi}\\}}\tilde{Q}_{j}^{i}(t)g_{\varepsilon}^{j}(\tilde{X}_{t}^{\xi},T-t)\textrm{d}t\right].$
(2.21)

By taking conditional expectation on the event that
$\tilde{X}_{T}^{\xi}=\eta$, we obtain that

$\displaystyle\mathbb{E}\left[1_{\\{T<T_{\xi}\\}}\tilde{Q}_{j}^{i}(T)W_{\varepsilon}^{j}(\tilde{X}_{T}^{\xi},0)\right]$
$\displaystyle=\int_{D}\mathbb{E}\left[\left.\tilde{Q}_{j}^{i}(T)1_{\\{T<T_{\xi}\\}}\right|\tilde{X}_{T}^{\xi}=\eta\right]W_{\varepsilon}^{j}(\eta,0)p_{-u_{T}}(0,\xi,T,\eta)\textrm{d}\eta$
$\displaystyle=\int_{D}\mathbb{E}\left[\left.\tilde{Q}_{j}^{i}(T)1_{\\{T<T_{\xi}\\}}\right|\tilde{X}_{T}^{\xi}=\eta\right]W_{\varepsilon}^{j}(\eta,0)p_{u}(0,\eta,T,\xi)\textrm{d}\eta$

and

$\displaystyle\mathbb{E}\left[\int_{0}^{T\wedge
T_{\xi}}\tilde{Q}_{j}^{i}(t)g_{\varepsilon}^{j}(\tilde{X}_{t}^{\xi},T-t)\textrm{d}t\right]$
$\displaystyle=\int_{0}^{T}\int_{D}\mathbb{E}\left[\left.\tilde{Q}_{j}^{i}(t)1_{\\{t<T_{\xi}\\}}g_{\varepsilon}^{j}(\tilde{X}_{t}^{\xi},T-t)\right|\tilde{X}_{T}^{\xi}=\eta\right]p_{-u_{T}}(0,\xi,T,\eta)\textrm{d}\eta\textrm{d}t$
$\displaystyle=\int_{0}^{T}\int_{D}\mathbb{E}\left[\left.\tilde{Q}_{j}^{i}(t)1_{\\{t<T_{\xi}\\}}g_{\varepsilon}^{j}(\tilde{X}_{t}^{\xi},T-t)\right|\tilde{X}_{T}^{\xi}=\eta\right]p_{u}(0,\eta,T,\xi)\textrm{d}\eta\textrm{d}t$

where $p_{-u_{T}}$ is the transition probability density function of the
diffusion $\tilde{X}^{\xi}$, which is a diffusion process with infinitesimal
generator $\nu\Delta-u_{T}\cdot\nabla$, where $u_{T}(x,t)=u(x,T-t)$. Since
$\nabla\cdot u=0$, we have that $p_{-u_{T}}(0,\xi,T,\eta)$ coincides with
$p_{u}(0,\eta,T,\xi)$. Therefore

$\displaystyle W_{\varepsilon}^{i}(\xi,T)$
$\displaystyle=\int_{D}\mathbb{E}\left[\left.\tilde{Q}_{j}^{i}(T)1_{\\{T<T_{\xi}\\}}\right|\tilde{X}_{T}^{\xi}=\eta\right]W_{\varepsilon}^{j}(\eta,0)p_{u}(0,\eta,T,\xi)\textrm{d}\eta$
$\displaystyle+\int_{0}^{T}\int_{D}\mathbb{E}\left[\left.\tilde{Q}_{j}^{i}(t)1_{\\{t<T_{\xi}\\}}g_{\varepsilon}^{j}(\tilde{X}_{t}^{\xi},T-t)\right|\tilde{X}_{T}^{\xi}=\eta\right]p_{u}(0,\eta,T,\xi)\textrm{d}\eta\textrm{d}t.$
(2.22)

Let us rewrite the representation (2.22) in terms of the distributions of the
Taylor diffusion. To this end, let us use $\tilde{\mathbb{P}}^{\xi}$ to denote
the law of $\tilde{X}^{\xi}$ and $\tilde{\mathbb{P}}^{\xi\rightarrow\eta}$ to
denote the conditional law, or the diffusion bridge measure of
$\tilde{X}^{\xi}$ given the terminal value that $\tilde{X}_{T}^{\xi}=\eta$.
The conditional law $\mathbb{\tilde{P}}^{\xi\rightarrow\eta}$ can be
considered as a probability measure on the path space
$C([0,T];\mathbb{R}^{3})$ canonically. For any continuous path
$\psi:[0,T]\rightarrow\mathbb{R}^{3}$, $\tilde{Q}(\psi;t)$ denotes the
solution to the following linear ordinary differential equations

$\frac{\textrm{d}}{\textrm{d}t}\tilde{Q}_{j}^{i}(\psi;t)=\tilde{Q}_{k}^{i}(\psi;t)q_{j}^{k}(\psi(t),T-t),\quad\tilde{Q}_{j}^{i}(\psi;0)=\delta_{j}^{i}$
(2.23)

for $i,j=1,2,3$. Then $\tilde{Q}_{j}^{i}(\tilde{X}^{\xi};t)$ gives rise a
version of the gauge functional defined above $\tilde{Q}_{j}^{i}(t)$.
Therefore, under the notations we just set up, the representation (2.22) may
be written as

$\displaystyle W_{\varepsilon}^{i}(\xi,T)$
$\displaystyle=\int_{D}\mathbb{\tilde{P}}^{\xi\rightarrow\eta}\left[\tilde{Q}_{j}^{i}(\psi;T)1_{\\{T<\zeta(\psi)\\}}\right]W_{\varepsilon}^{j}(\eta,0)p_{u}(0,\eta,T,\xi)\textrm{d}\eta$
$\displaystyle+\int_{0}^{T}\int_{D}\mathbb{\tilde{P}}^{\xi\rightarrow\eta}\left[\tilde{Q}_{j}^{i}(\psi;t)1_{\\{t<\zeta(\psi)\\}}g_{\varepsilon}^{j}(\psi(t),T-t)\right]p_{u}(0,\eta,T,\xi)\textrm{d}\eta\textrm{d}t$
(2.24)

where

$\zeta(\psi)=\inf\left\\{t\geq 0:\psi(t)\notin D\right\\}.$

Next we apply the fundamental duality for conditional laws: since $\nabla\cdot
u=0$ on $\mathbb{R}^{3}$ in the distribution sense, so that

$\mathbb{\tilde{P}}^{\xi\rightarrow\eta}=\mathbb{P}^{\eta\rightarrow\xi}\circ\tau_{T},$

cf. [37]. Let $Q(\psi,T;t)=\tilde{Q}(\psi\circ\tau_{T};T-t)$. Since

$\tilde{Q}_{j}^{i}(\psi;s)=\delta_{ij}+\int_{0}^{s}\tilde{Q}_{k}^{i}(\psi;r)q_{j}^{k}(\psi(r),T-r)\textrm{d}r.$
(2.25)

Apply this to $\psi\circ\tau_{T}$ one obtains that

$\displaystyle\tilde{Q}_{j}^{i}(\psi\circ\tau_{T};s)$
$\displaystyle=\delta_{ij}+\int_{0}^{s}\tilde{Q}_{k}^{i}(\psi\circ\tau_{T};r)q_{j}^{k}(\psi(T-r),T-r)\textrm{d}r$
$\displaystyle=\delta_{ij}+\int_{T-s}^{T}\tilde{Q}_{k}^{i}(\psi\circ\tau_{T};T-r)q_{j}^{k}(\psi(r),r)\textrm{d}r$

so that

$\tilde{Q}_{j}^{i}(\psi\circ\tau_{T};T-s)=\delta_{ij}+\int_{s}^{T}\tilde{Q}_{k}^{i}(\psi\circ\tau_{T};T-r)q_{j}^{k}(\psi(r),r)\textrm{d}r.$

This means that $t\mapsto Q(\psi,T;t)$ is the unique solution to the ordinary
differential equation:

$dQ_{j}^{i}(\psi,T;t)=-Q_{k}^{i}(\psi,T;t)q_{j}^{k}(\psi(t),t)dt,\quad
Q_{j}^{i}(\psi,T;T)=\delta_{ij}.$

By definition, $\tilde{Q}(\psi;t)=Q(\psi\circ\tau_{T},T;T-t)$ for every
$t\in[0,T].$ Hence, by the duality of the conditional laws, we may rewrite
(2.24) in terms of the law of the Taylor diffusion:

$\displaystyle W_{\varepsilon}^{i}(\xi,T)$
$\displaystyle=\int_{D}\mathbb{P}^{\eta\rightarrow\xi}\left[Q_{j}^{i}(\psi,T;0)1_{\\{T<\zeta(\psi\circ\tau_{T})\\}}\right]W_{\varepsilon}^{j}(\eta,0)p_{u}(0,\eta,T,\xi)\textrm{d}\eta$
$\displaystyle+\int_{0}^{T}\int_{D}\mathbb{P}^{\eta\rightarrow\xi}\left[Q_{j}^{i}(\psi,T;T-t)1_{\\{t<\zeta(\psi\circ\tau_{T})\\}}g_{\varepsilon}^{j}(\psi(T-t),T-t)\right]p_{u}(0,\eta,T,\xi)\textrm{d}\eta\textrm{d}t$
$\displaystyle=\int_{D}\mathbb{P}^{\eta\rightarrow\xi}\left[Q_{j}^{i}(\psi,T;0)1_{\\{T<\zeta(\psi\circ\tau_{T})\\}}\right]W_{\varepsilon}^{j}(\eta,0)p_{u}(0,\eta,T,\xi)\textrm{d}\eta$
$\displaystyle+\int_{0}^{T}\int_{D}\mathbb{P}^{\eta\rightarrow\xi}\left[Q_{j}^{i}(\psi,T;t)1_{\\{T-t<\zeta(\psi\circ\tau_{T})\\}}g_{\varepsilon}^{j}(\psi(t),t)\right]p_{u}(0,\eta,T,\xi)\textrm{d}\eta\textrm{d}t$

which yields the claim. ∎

We are now in a position to derive an important stochastic representation
formula.

######  Theorem 2.2.

Let $X^{\eta}$ be the Taylor diffusion:

$\textrm{d}X_{t}^{\eta}=u(X_{t}^{\eta},t)\textrm{d}t+\sqrt{2\nu}\textrm{d}B_{t},\quad
X_{0}^{\eta}=\eta$ (2.26)

for each $\eta\in\mathbb{R}^{3}$. The following stochastic representation
holds:

$\displaystyle u^{i}(x,t)$
$\displaystyle=\varepsilon^{ijk}\int_{D}K^{j}(x,\eta)\sigma_{\varepsilon}^{k}(\eta,t)\textrm{d}\eta$
$\displaystyle+\varepsilon^{ijk}\int_{D}\mathbb{E}\left[Q_{l}^{k}(\eta,t;0)1_{\\{t<\zeta(X^{\eta}\circ\tau_{t})\\}}K^{j}(x,X_{t}^{\eta})\right]W_{\varepsilon}^{l}(\eta,0)\textrm{d}\eta$
$\displaystyle+\varepsilon^{ijk}\int_{0}^{t}\int_{D}\mathbb{E}\left[Q_{l}^{k}(\eta,t;s)1_{\\{t-s<\zeta(X^{\eta}\circ\tau_{t})\\}}K^{j}(x,X_{t}^{\eta})g_{\varepsilon}^{l}(X_{s}^{\eta},s)\right]\textrm{d}\eta\textrm{d}s$

for any $t>0$ and $x\in D$, where $i=1,2,3$.

###### Proof.

Recall the Biot-Savart law

$u(x,t)=\int_{D}K(x,\xi)\wedge\omega(\xi,t)\textrm{d}\xi$

where $K$ is the singular kernel of the half-space $x_{3}>0$. By Theorem 2.1

$\displaystyle\omega^{k}(\xi,t)$
$\displaystyle=\sigma_{\varepsilon}^{k}(\xi,t)+\int_{D}\mathbb{P}^{\eta\rightarrow\xi}\left[Q_{l}^{k}(\eta,t;0)1_{\\{t<\zeta(X^{\eta}\circ\tau_{t})\\}}\right]W_{\varepsilon}^{l}(\eta,0)p_{u}(0,\eta,t,\xi)\textrm{d}\eta$
$\displaystyle+\int_{0}^{t}\int_{D}\mathbb{P}^{\eta\rightarrow\xi}\left[Q_{l}^{k}(\eta,t;s)1_{\\{t-s<\zeta(X^{\eta}\circ\tau_{t})\\}}g_{\varepsilon}^{l}(X_{s}^{\eta},s)\right]p_{u}(0,\eta,t,\xi)\textrm{d}\eta\textrm{d}s.$
(2.27)

Integrating with $\varepsilon^{ijk}K^{j}(x,\cdot)$ and using Fubini theorem,
we obtain

$\displaystyle u^{i}(x,t)$
$\displaystyle=\varepsilon^{ijk}\int_{D}K^{j}(x,\eta)\sigma_{\varepsilon}^{k}(\eta,t)\textrm{d}\eta$
$\displaystyle+\varepsilon^{ijk}\int_{D}\mathbb{E}\left[K^{j}(x,X_{t}^{\eta})Q_{l}^{k}(\eta,t;0)1_{\\{t<\zeta(X^{\eta}\circ\tau_{t})\\}}\right]W_{\varepsilon}^{l}(\eta,0)\textrm{d}\eta$
$\displaystyle+\varepsilon^{ijk}\int_{0}^{t}\int_{D}\mathbb{E}\left[Q_{l}^{k}(\eta,t;s)1_{\\{t-s<\zeta(X^{\eta}\circ\tau_{t})\\}}K^{j}(x,X_{t}^{\eta})g_{\varepsilon}^{l}(X_{s}^{\eta},s)\right]\textrm{d}\eta\textrm{d}s$

which completes the proof. ∎

###  2.2 Two dimensional wall-bounded flows

There is a simplified version for two dimensional case. For two dimensional
flows, the non-linear stretching term $A^{i}_{j}\omega^{j}$ vanishes
identically, so that we can take $Q^{i}_{j}=\delta_{ij}$. Therefore we have
the following two dimensional representation formula.

######  Theorem 2.3.

Suppose $d=2$, so that $\omega=\frac{\partial}{\partial
x_{1}}u^{2}-\frac{\partial}{\partial x_{2}}u^{1}$ is a scalar function. Then

$\displaystyle\omega(\xi,t)$
$\displaystyle=\sigma_{\varepsilon}(\xi,t)+\int_{D}\mathbb{P}^{\eta\rightarrow\xi}\left[t<\zeta(X^{\eta}\circ\tau_{t})\right]W_{\varepsilon}(\eta,0)p_{u}(0,\eta,t,\xi)\textrm{d}\eta$
$\displaystyle+\int_{0}^{t}\int_{D}\mathbb{E}^{\eta\rightarrow\xi}\left[1_{\\{t-s<\zeta(X^{\eta}\circ\tau_{t})\\}}g_{\varepsilon}(X_{s}^{\eta},s)\right]p_{u}(0,\eta,t,\xi)\textrm{d}\eta\textrm{d}s$
(2.28)

for every $\xi\in D$ and $t>0$, where $\zeta(\psi)=\inf\left\\{s:\psi(s)\notin
D\right\\}$ and $X^{\eta}$ is the Taylor diffusion:

$\textrm{d}X_{t}^{\eta}=u(X_{t}^{\eta},t)\textrm{d}t+\sqrt{2\nu}\textrm{d}B_{t},\quad
X_{0}^{\eta}=\eta$ (2.29)

for every $\eta\in\mathbb{R}^{2}$.

As a consequence, we have the following representation for the velocity $u$ in
two dimensional case.

######  Theorem 2.4.

Suppose $d=2$, and let $K=(K^{1},K^{2})$ be the Biot-Savart kernel for the
half-plane $D:x_{2}>0$. Then

$\displaystyle u^{i}(x,t)$
$\displaystyle=\int_{D}K^{i}(x,\eta)\sigma_{\varepsilon}(\eta,t)\textrm{d}\eta+\int_{D}\mathbb{E}\left[K^{i}(x,X_{t}^{\eta})1_{\left\\{t<\zeta(X^{\eta}\circ\tau_{t})\right\\}}\right]W_{\varepsilon}(\eta,0)\textrm{d}\eta$
$\displaystyle+\int_{0}^{t}\int_{D}\mathbb{E}\left[1_{\\{t-s<\zeta(X^{\eta}\circ\tau_{t})\\}}K^{i}(x,X_{t}^{\eta})g_{\varepsilon}(X_{s}^{\eta},s)\right]\textrm{d}\eta\textrm{d}s$
(2.30)

for every $x\in D$ and $t>0$, where $\zeta(\psi)=\inf\left\\{s:\psi(s)\notin
D\right\\}$ and $X^{\eta}$ is the Taylor diffusion:

$\textrm{d}X_{t}^{\eta}=u(X_{t}^{\eta},t)\textrm{d}t+\sqrt{2\nu}\textrm{d}B_{t},\quad
X_{0}^{\eta}=\eta$ (2.31)

for every $\eta\in\mathbb{R}^{2}$.

##  3 2D wall-bounded flows

In the subsequent parts we aim to present a stochastic formulation of a
viscous fluid flow, with kinematic viscosity constant $\nu>0$, which passes
through obstacles, such as a thin plane plate and a wedge. These flows have
been studied as important and classical examples which demonstrate boundary
layer phenomena. Turbulence may be build up near the solid obstacle when the
Reynolds number becomes large.

For simplicity we take a two dimensional model, although we admit that the
three dimensional model is more sophisticated and will be studied in a future
work. Therefore the fluid flow in question is described by its velocity, a
time dependent vector field $u=(u^{1},u^{2})$ in the domain $D$.

###  3.1 A general case

Let $D\subset\mathbb{R}^{2}$ be a general simply connected domain which is not
the whole plane. For such domain, according to Riemann mapping theorem, there
is a one-to-one and onto conformal mapping $T:D\rightarrow\mathbb{H}$ where
$\mathbb{H}$ is the upper half-plane (the mapping is unique up to a rotation
at a point whose image is assigned to $i$ for example). Then the Green
function for such domain $D$ is given by

$G_{D}(x,y)=\Gamma(T(x),T(y))-\Gamma(T(x),\overline{T(y)})$ (3.1)

for $x,y\in D$, where $\Gamma(x,y)=\frac{1}{2\pi}\log|x-y|$. While the
boundary $\partial D$ may be complicated, so some care is needed, although
formally $G$ satisfies the Dirichlet boundary condition.

Notice also that $\omega=\nabla\wedge u$ and $\nabla\cdot u=0$ imply

$\Delta u^{1}=-\frac{\partial}{\partial x_{2}}\omega,\quad\Delta
u^{2}=\frac{\partial}{\partial x_{1}}\omega\quad\text{ in }D,$ (3.2)

and, as $u^{i}$ are subject to the no-slip condition, the Green formula
implies the Biot-Savart law

$u^{i}(x,t)=\int_{D}K^{i}(x,y)\omega(y,t)\textrm{d}y,$ (3.3)

where the kernel $K$ is given by

$K^{1}(x,y)=\frac{\partial}{\partial y_{2}}G_{D}(x,y),\quad
K^{2}(x,y)=-\frac{\partial}{\partial y_{1}}G_{D}(x,y).$ (3.4)

Hence, given the Green function as in (3.1), we might state the following
lemma.

######  Lemma 3.1.

The kernel $K$ can be computed as follows

$\begin{split}K^{1}(x,y)&=k^{-}(x,y)\cdot\frac{\partial T}{\partial
y_{2}}(y)-k^{+}(x,y)\cdot\frac{\partial T}{\partial y_{2}}(y),\\\
K^{2}(x,y)&=-k^{-}(x,y)\cdot\frac{\partial T}{\partial
y_{1}}(y)+k^{+}(x,y)\cdot\frac{\partial T}{\partial y_{1}}(y),\end{split}$
(3.5)

where the vector functions $k^{-}(x,y),k^{+}(x,y)$ are given by

$\begin{split}k^{-}(x,y)=\frac{1}{2\pi}\frac{(T^{1}(y)-T^{1}(x),T^{2}(y)-T^{2}(x))}{(T^{1}(y)-T^{1}(x))^{2}+(T^{2}(y)-T^{2}(x))^{2}},\\\
k^{+}(x,y)=\frac{1}{2\pi}\frac{(T^{1}(y)-T^{1}(x),T^{2}(y)+T^{2}(x))}{(T^{1}(y)-T^{1}(x))^{2}+(T^{2}(y)+T^{2}(x))^{2}},\end{split}$
(3.6)

and $\frac{\partial T}{\partial y_{1}}(y)$ and $\frac{\partial T}{\partial
y_{2}}(y)$ denote the vectors $(\frac{\partial T^{1}}{\partial
y_{1}}(y),\frac{\partial T^{2}}{\partial y_{1}}(y))$ and $(\frac{\partial
T^{1}}{\partial y_{2}}(y),\frac{\partial T^{2}}{\partial y_{2}}(y))$
respectively.

As for the numerical simulations we also need to compute the boundary stress
$\theta$, we do so by differentiating the representation (3.3) with respect to
$x$ (and taking limits as $x$ approaches the boundary). Note that in practice
we have to desingularise the kernels (see Section 5 for details), which
insures we can differentiate the singular integral representation. Therefore,
we can write the representation, according to the lemma above, in terms of
derivatives of $k^{-},k^{+}$ given explicitly in the following.

######  Lemma 3.2.

The derivatives $\dfrac{\partial k^{-}}{\partial x},\dfrac{\partial
k^{+}}{\partial x}$ are given by

$\begin{split}\dfrac{\partial(k^{-})^{1}}{\partial
x_{i}}(x,y)=\frac{1}{2\pi}\Bigg{(}\frac{-\frac{\partial T^{1}}{\partial
x_{i}}(x)}{(T^{1}(y)-T^{1}(x))^{2}+(T^{2}(y)-T^{2}(x))^{2}}+2(T^{1}(y)-T^{1}(x))\\\
\cdot\frac{(T^{1}(y)-T^{1}(x))\frac{\partial T^{1}}{\partial
x_{i}}(x)+(T^{2}(y)-T^{2}(x))\frac{\partial T^{2}}{\partial
x_{i}}(x)}{\big{(}(T^{1}(y)-T^{1}(x))^{2}+(T^{2}(y)-T^{2}(x))^{2}\big{)}^{2}}\Bigg{)},\\\
\dfrac{\partial(k^{-})^{2}}{\partial
x_{i}}(x,y)=\frac{1}{2\pi}\Bigg{(}\frac{-\frac{\partial T^{2}}{\partial
x_{i}}(x)}{(T^{1}(y)-T^{1}(x))^{2}+(T^{2}(y)-T^{2}(x))^{2}}+2(T^{2}(y)-T^{2}(x))\\\
\cdot\frac{(T^{1}(y)-T^{1}(x))\frac{\partial T^{1}}{\partial
x_{i}}(x)+(T^{2}(y)-T^{2}(x))\frac{\partial T^{2}}{\partial
x_{i}}(x)}{\big{(}(T^{1}(y)-T^{1}(x))^{2}+(T^{2}(y)-T^{2}(x))^{2}\big{)}^{2}}\Bigg{)},\end{split}$
(3.7)

and

$\begin{split}\dfrac{\partial(k^{+})^{1}}{\partial
x_{i}}(x,y)=\frac{1}{2\pi}\Bigg{(}\frac{-\frac{\partial T^{1}}{\partial
x_{i}}(x)}{(T^{1}(y)-T^{1}(x))^{2}+(T^{2}(y)+T^{2}(x))^{2}}+2(T^{1}(y)-T^{1}(x))\\\
\cdot\frac{(T^{1}(y)-T^{1}(x))\frac{\partial T^{1}}{\partial
x_{i}}(x)-(T^{2}(y)+T^{2}(x))\frac{\partial T^{2}}{\partial
x_{i}}(x)}{\big{(}(T^{1}(y)-T^{1}(x))^{2}+(T^{2}(y)+T^{2}(x))^{2}\big{)}^{2}}\Bigg{)},\\\
\dfrac{\partial(k^{+})^{2}}{\partial
x_{i}}(x,y)=\frac{1}{2\pi}\Bigg{(}\frac{\frac{\partial T^{2}}{\partial
x_{i}}(x)}{(T^{1}(y)-T^{1}(x))^{2}+(T^{2}(y)+T^{2}(x))^{2}}+2(T^{2}(y)+T^{2}(x))\\\
\cdot\frac{(T^{1}(y)-T^{1}(x))\frac{\partial T^{1}}{\partial
x_{i}}(x)-(T^{2}(y)+T^{2}(x))\frac{\partial T^{2}}{\partial
x_{i}}(x)}{\big{(}(T^{1}(y)-T^{1}(x))^{2}+(T^{2}(y)+T^{2}(x))^{2}\big{)}^{2}}\Bigg{)},\end{split}$
(3.8)

where $(k^{\pm})^{i}$ denote the components of $k^{\pm}$, i.e.
$k^{\pm}=((k^{\pm})^{1},(k^{\pm})^{2})$.

###  3.2 Flows past a thin plate

We consider a flow in the domain
$D=\mathbb{R}^{2}\setminus\left\\{(x_{1},0):x_{1}\geq 0\right\\}$, where the
boundary $\partial D$: $x_{2}=0$ and $x_{1}\geq 0$ is an obstacle in the flow.
The fluid flow has to satisfy the no-slip condition

$\lim_{x_{2}\uparrow 0}u((x_{1},x_{2}),t)=\lim_{x_{2}\downarrow
0}u((x_{1},x_{2}),t)=0,\quad\text{when }x_{1}\geq 0$ (3.9)

for $t>0$. Therefore automatically $u(x,t)$ is extended to the whole plane by
setting $u(x,t)=0$ for $x\in\partial D$. This extension has the following nice
property: if $\nabla\cdot u=0$ in $D$, then $\nabla\cdot u=0$ in
$\mathbb{R}^{2}$ in the sense of distributions.

To find the Green function we use the conformal principle. The domain $D$ is
conformally equivalent to $\mathbb{H}$, it is realised by the conformal
transform $T^{-1}:x\rightarrow x^{2}$ (consider $x=(x_{1},x_{2})=x_{1}+ix_{2}$
as a complex coordinate) from $\mathbb{H}$ one-to-one and onto $D$. In terms
of real coordinates

$T^{-1}:(x_{1},x_{2})\rightarrow(x_{1}^{2}-x_{2}^{2},2x_{1}x_{2})$

and its inverse

$T:(x_{1},x_{2})\rightarrow\left(\textrm{sgn}(x_{2})\sqrt{\frac{1}{2}\left(|x|+x_{1}\right)},\sqrt{\frac{1}{2}\left(|x|-x_{1}\right)}\right)$
(3.10)

where $|x|=\sqrt{x_{1}^{2}+x_{2}^{2}}$.

Recall that the Green function for $D$ is given by

$G(x,y)=\Gamma(Tx,Ty)-\Gamma(Tx,\overline{Ty}),$

which vanishes for $x$ such that

$Tx=\overline{Tx},$

which is equivalent to that $|x|-x_{1}=0$, i.e. $x_{1}\geq 0$ and $x_{2}=0$.
To compute the kernel $K$ in this case, we explicitly find the partial
derivatives of the map $T$ and use Lemma 3.1.

######  Lemma 3.3.

Define $T=T^{1}+iT^{2}$ by (3.10). Its partial derivatives are given by

$\begin{split}\frac{\partial T^{1}}{\partial
x_{1}}(x)&=\textrm{sgn}(x_{2})\frac{\sqrt{\frac{1}{2}\left(|x|+x_{1}\right)}}{2|x|},\\\
\frac{\partial T^{2}}{\partial
x_{1}}(x)&=-\frac{\sqrt{\frac{1}{2}\left(|x|-x_{1}\right)}}{2|x|},\end{split}$
(3.11)

and

$\begin{split}\frac{\partial T^{1}}{\partial
x_{2}}(x)&=\textrm{sgn}(x_{2})\frac{x_{2}}{4|x|\sqrt{\frac{1}{2}\left(|x|+x_{1}\right)}},\\\
\frac{\partial T^{2}}{\partial
x_{2}}(x)&=\frac{x_{2}}{4|x|\sqrt{\frac{1}{2}\left(|x|-x_{1}\right)}}.\end{split}$
(3.12)

###  3.3 Flows past a wedge obstacle

Suppose the wedge with angle $2\alpha$ at its tip (where
$\alpha\in(0,\frac{\pi}{2})$) is modelled by

$\Lambda=\left\\{x=r\textrm{e}^{\textrm{i}\theta}:r\geq 0\textrm{ and
}-\alpha\leq\theta\leq\alpha\right\\}$

so that the fluid occupies the region $D=\mathbb{R}^{2}\setminus\Lambda$. It
is again conformally equivalent to the upper half space $\mathbb{H}$, the
conformal mapping $T:D\rightarrow\mathbb{H}$ is given explicitly as follows:

$T(x)=\left(\textrm{e}^{-\textrm{i}\alpha}x\right)^{\beta}=\exp\left[\beta\ln|x|+\textrm{i}\beta\arg(\textrm{e}^{-\textrm{i}\alpha}x)\right]\quad\textrm{
for }x\in D$

where $0<\arg x<2\pi$,

$\beta=\frac{\pi}{2(\pi-\alpha)}\quad\textrm{ and
}\quad\alpha\beta=\frac{\pi\alpha}{2(\pi-\alpha)}$

and

$x^{\beta}=|x|^{\beta}\exp\left(\textrm{i}\beta\arg x\right)\quad\textrm{ with
}0<\arg x<2\pi.$

Therefore

$T(x)=\exp\left(\textrm{i}\beta\arg\left(e^{-\textrm{i}\alpha}x\right)\right)|x|^{\beta}\quad\textrm{
for }x\in D.$

In terms of real coordinates we have

$T(x)=\left(|x|^{\beta}\cos\left(\beta\arg\left(e^{-\textrm{i}\alpha}x\right)\right),|x|^{\beta}\sin\left(\beta\arg\left(e^{-\textrm{i}\alpha}x\right)\right)\right)$

for $x=(x_{1},x_{2})\in D$, where $\arg x\in(0,2\pi)$, so that

$\arg(x_{1},x_{2})=\begin{cases}\arctan\frac{x_{2}}{x_{1}}&\textrm{if
}x_{1}>0,x_{2}\geq 0,\\\ \frac{\pi}{2}&\textrm{if }x_{1}=0,x_{2}>0,\\\
\pi+\arctan\frac{x_{2}}{x_{1}}&\textrm{if }x_{1}<0,\\\
\frac{3\pi}{2}&\textrm{if }x_{1}=0,x_{2}<0,\\\
2\pi+\arctan\frac{x_{2}}{x_{1}}&\textrm{if }x_{1}>0,x_{2}\leq 0.\end{cases}$

Since $\alpha\in(0,\frac{\pi}{2})$ and $\arg x\in(\alpha,2\pi-\alpha)$ for
$x\in D$, it follows that

$\arg\left(e^{-\textrm{i}\alpha}x\right)=\arg x-\alpha,$

hence the real and imaginary parts $T=T^{1}+iT^{2}$ can be written as

$T^{1}(x)=|x|^{\beta}\cos{(\beta(\arg x-\alpha))},\quad
T^{2}(x)=|x|^{\beta}\sin{(\beta(\arg x-\alpha))}.$ (3.13)

Therefore the Green function for this domain is given by

$G(x,y)=\Gamma(Tx,Ty)-\Gamma(Tx,\overline{Ty})$

which vanishes for $x$ such that

$Tx=\overline{Tx}$

that is

$|x|^{\beta}\sin\left(\beta\arg\left(e^{-\textrm{i}\alpha}x\right)\right)=-|x|^{\beta}\sin\left(\beta\arg\left(e^{-\textrm{i}\alpha}x\right)\right)$

so that either $x=(0,0)$ or

$\sin\left(\beta\arg(e^{-\textrm{i}\alpha}x)\right)=0.$

That is

$\beta\arg(e^{-\textrm{i}\alpha}x)=n\pi$

so that

$\arg(e^{-\textrm{i}\alpha}x)=2n(\pi-\alpha).$

That is

$\arg x-\alpha=2n(\pi-\alpha)\quad\textrm{ and }0<\arg x-\alpha<2\pi.$

The requirement

$0<2n(\pi-\alpha)<2\pi$

gives rise to the constraint that $n=0,1$ only that is, $\arg x=\alpha$ and
$\arg x=2\pi-\alpha$ which are the edges of the wedge.

For computing the Biot-Savart kernel $K$ using Lemma 3.1, we write the partial
derivatives of the map $T$.

######  Lemma 3.4.

Let $T=T^{1}+iT^{2}$ be as in (3.13), then its partial derivatives are given
by

$\begin{split}\frac{\partial T^{1}}{\partial
x_{1}}(x)=\beta|x|^{\beta-2}(x_{1}\cos{(\beta(\arg
x-\alpha))}+x_{2}\sin{(\beta(\arg x-\alpha))}),\\\ \frac{\partial
T^{2}}{\partial x_{1}}(x)=\beta|x|^{\beta-2}(x_{1}\sin{(\beta(\arg
x-\alpha))}-x_{2}\cos{(\beta(\arg x-\alpha))}),\end{split}$ (3.14)

and

$\begin{split}\frac{\partial T^{1}}{\partial
x_{2}}(x)=\beta|x|^{\beta-2}(x_{2}\cos{(\beta(\arg
x-\alpha))}-x_{1}\sin{(\beta(\arg x-\alpha))}),\\\ \frac{\partial
T^{2}}{\partial x_{2}}(x)=\beta|x|^{\beta-2}(x_{2}\sin{(\beta(\arg
x-\alpha))}+x_{1}\cos{(\beta(\arg x-\alpha))}).\end{split}$ (3.15)

##  4 Limiting stochastic representations

The goal of this section is to derive representations similar to those in
Section 2 for different domains described above. Note that so far we have a
stochastic representation for the velocity $u$, e.g., the one given in Theorem
2.4, for the half-plane. This representation contains the time-integral with
the term $g_{\varepsilon}$ given by

$g_{\varepsilon}=G-\frac{\partial\sigma_{\varepsilon}}{\partial
t}-(u\cdot\nabla)\sigma_{\varepsilon}+\nu\Delta\sigma_{\varepsilon},$ (4.1)

where $\varepsilon$ defines the width of the cutoff modification in
$\sigma_{\varepsilon}$. It is therefore interesting theoretically to consider
the limits of the representations as the cutoff width parameter $\varepsilon$
converges to zero. It turns out that many terms in $g_{\varepsilon}$ do not
contribute to the limit, we are able to simplify the expressions used for
simulations in Section 5.

###  4.1 Flat-plate case

Let us deal with the flows passing a flat-plate first, that is
$D=\\{x=(x_{1},x_{2}):x_{1}<0\textrm{ or }x_{2}\neq 0\\}$. The flat-plate is
modelled by the boundary $\partial D=\\{x_{2}=0,x_{1}\geq 0\\}$, and the flow
is split into upper flow and lower flow after hitting the plate. The velocity
$u$ has to satisfy the no-slip condition, that is,

$\lim_{x_{2}\rightarrow 0\pm}u(x_{1},x_{2},t)=0\quad\textrm{ for }x_{1}\geq
0\textrm{ and }t>0.$ (4.2)

Therefore there is non-trivial stress on both sides of the plate. The stress
at the upper side of the plate is denoted by $\theta_{+}$ and the one at the
lower side by $\theta_{-}$, that is

$\theta_{\pm}(x_{1},t)=\begin{cases}-\left.\frac{\partial}{\partial
x_{2}}u^{1}(x_{1},x_{2},t)\right|_{x_{2}=0\pm}&\textrm{ for }x_{1}\geq 0,\\\
0&\textrm{ for }x_{1}<0.\end{cases}$ (4.3)

We make a technical assumption that $\theta_{+}$ and $\theta_{-}$ have
continuous derivatives up to second order except at $x_{1}=0$. We also assume
that the limit of $-\frac{\partial}{\partial x_{2}}u^{1}(x_{1},x_{2},t)$
exists as $(x_{1},x_{2})\to(0,0)$ in $D$, therefore
$\theta_{+}(0,t)=\theta_{-}(0,t)$, for all $t\geq 0$, denoted $\theta(t)$ in
this case. Let $\phi:\mathbb{R}\rightarrow[0,1]$ be a cut-off function which
is smooth on $(0,\infty)$ such that $\phi(r)=0$ for $r<0$, $\phi(r)=1$ for
$0\leq r<1/3$ and $\phi(r)=0$ for $r\geq 2/3$. Define

$\sigma_{\varepsilon}(x_{1},x_{2},t)=\theta_{+}(x_{1},t)\phi(x_{2}/\varepsilon)+\theta_{-}(x_{1},t)\phi(-x_{2}/\varepsilon)+\theta(t)\phi(|x|/\varepsilon)1_{\\{x_{1}<0\\}}$
(4.4)

for $x=(x_{1},x_{2})\in\mathbb{R}^{2}$.

Let $W_{\varepsilon}=\omega-\sigma_{\varepsilon}$. Then

$\left(\frac{\partial}{\partial
t}+u\cdot\nabla-\nu\Delta\right)W_{\varepsilon}=G+\rho_{\varepsilon}\quad\textrm{
in }D,$ (4.5)

satisfying the boundary condition $\left.W_{\varepsilon}\right|_{\partial
D}=0$ in the sense that

$\lim_{x_{2}\rightarrow 0\pm}W_{\varepsilon}(x_{1},x_{2},t)=0\quad\textrm{ for
}x_{1}\geq 0,$ (4.6)

where

$\rho_{\varepsilon}=\nu\Delta\sigma_{\varepsilon}-\frac{\partial\sigma_{\varepsilon}}{\partial
t}-(u\cdot\nabla)\sigma_{\varepsilon}.$

Taking into account that $\sigma_{\varepsilon}$ is given as in (4.4), we
compute

$\displaystyle\rho_{\varepsilon}(x,t)$
$\displaystyle=\frac{\nu}{\varepsilon^{2}}\phi^{\prime\prime}(x_{2}/\varepsilon)\theta_{+}(x_{1},t)-\frac{1}{\varepsilon}\phi^{\prime}(x_{2}/\varepsilon)u^{2}(x,t)\theta_{+}(x_{1},t)$
$\displaystyle+\phi(x_{2}/\varepsilon)\left(\nu\frac{\partial^{2}\theta_{+}}{\partial
x_{1}^{2}}(x_{1},t)-\frac{\partial\theta_{+}}{\partial
t}(x_{1},t)\right)-\phi(x_{2}/\varepsilon)u^{1}(x,t)\frac{\partial\theta_{+}}{\partial
x_{1}}(x_{1},t)$
$\displaystyle+\frac{\nu}{\varepsilon^{2}}\phi^{\prime\prime}(-x_{2}/\varepsilon)\theta_{-}(x_{1},t)+\frac{1}{\varepsilon}\phi^{\prime}(-x_{2}/\varepsilon)u^{2}(x,t)\theta_{-}(x_{1},t)$
$\displaystyle+\phi(-x_{2}/\varepsilon)\left(\nu\frac{\partial^{2}\theta_{-}}{\partial
x_{1}^{2}}(x_{1},t)-\frac{\partial\theta_{-}}{\partial
t}(x_{1},t)\right)-\phi(-x_{2}/\varepsilon)u^{1}(x,t)\frac{\partial\theta_{-}}{\partial
x_{1}}(x_{1},t)$
$\displaystyle+\Big{(}-\frac{1}{\varepsilon}\theta(t)\phi^{\prime}(|x|/\varepsilon)\frac{u^{1}(x_{1},x_{2})x_{1}+u^{2}(x_{1},x_{2})x_{2}}{|x|}$
$\displaystyle-\theta^{\prime}(t)\phi(|x|/\varepsilon)+\frac{\nu}{\varepsilon^{2}}\theta(t)\phi^{\prime\prime}(|x|/\varepsilon)+\frac{\nu}{\varepsilon}\theta(t)\phi^{\prime}(|x|/\varepsilon)\frac{1}{|x|}\Big{)}1_{\\{x_{1}<0\\}}$
(4.7)

for any $x\in D$. The initial data for $W_{\varepsilon}$ is given as follows

$W_{\varepsilon}(x,0)=\omega_{0}(x_{1},x_{2})-\theta_{+}(x_{1},0)\phi(x_{2}/\varepsilon)-\theta_{-}(x_{1},0)\phi(-x_{2}/\varepsilon)-\theta(0)\phi(|x|/\varepsilon)1_{\\{x_{1}<0\\}},$
(4.8)

for $x\in D$.

The velocity $u(x,t)$ is extended to the whole space $\mathbb{R}^{2}$
trivially by defining $u(x,t)=0$ if $x\notin D$. Then $\nabla\cdot u=0$ on
$\mathbb{R}^{2}$ in distribution sense. Let $p(s,\xi,t,\eta)$ denote the
transition probability density function of the diffusion process with
infinitesimal generator $\nu\Delta+u\cdot\nabla$, and $p^{D}(s,\xi,t,\eta)$ be
the transition probability density of the same diffusion killed on hitting the
boundary $\partial D$. Since $u$ is divergence-free, $p^{D}(s,\xi,t,\eta)$
coincides with the Green function of the Dirichlet problem associated with the
heat operator $\nu\Delta-u\cdot\nabla-\frac{\partial}{\partial t}$ on $D$,
therefore (cf. [20, Chapter 1, Theorem 12])

$\omega(y,t)=\sigma_{\varepsilon}(y,t)+\int_{D}p^{D}(0,\xi,t,y)W_{\varepsilon}(\xi,0)\textrm{d}\xi+\int_{0}^{t}\int_{D}p^{D}(s,\xi,t,y)(G(\xi,s)+\rho_{\varepsilon}(\xi,s))\textrm{d}\xi\textrm{d}s$
(4.9)

for $y\in D$ and $t>0$.

######  Lemma 4.1.

The following integral representation holds

$\displaystyle\omega(y,t)$
$\displaystyle=(\omega_{0}(y_{1},0+,t)+\omega_{0}(y_{1},0-,t))1_{\\{y\in\partial
D\\}}+\int_{D}p^{D}(0,\xi,t,y)\omega_{0}(\xi)\textrm{d}\xi$
$\displaystyle+\int_{0}^{t}\int_{D}p^{D}(s,\xi,t,y)G(\xi,s)\textrm{d}\xi\textrm{d}s$
$\displaystyle+\nu\int_{0}^{t}\int_{0}^{\infty}\frac{\partial}{\partial\xi_{2}}p^{D}(s,(\xi_{1},0+),t,y)\theta_{+}(\xi_{1},s)\textrm{d}\xi_{1}\textrm{d}s$
$\displaystyle-\nu\int_{0}^{t}\int_{0}^{\infty}\frac{\partial}{\partial\xi_{2}}p^{D}(s,(\xi_{1},0-),t,y)\theta_{-}(\xi_{1},s)\textrm{d}\xi_{1}\textrm{d}s$
(4.10)

for $y\in D$ and $t>0$.

###### Proof.

Following [38], we obtain the result by letting $\varepsilon\downarrow 0$ in
(4.9) — one can also take the same cutoff function $\phi$ given by

$\phi(r)=\begin{cases}1&\text{ if }0\leq r\leq\frac{1}{3},\\\
\frac{1}{2}+54\Big{(}r-\frac{1}{2}\Big{)}^{3}-\frac{9}{2}\Big{(}r-\frac{1}{2}\Big{)}&\text{
if }\frac{1}{3}<r\leq\frac{2}{3},\\\ 0&\text{ if }r>\frac{2}{3}.\end{cases}$
(4.11)

Taking (4.4) to the limit, we have

$\sigma_{\varepsilon}(y,t)\to\theta_{+}(y_{1},t)1_{\\{y_{2}=0\\}}+\theta_{-}(y_{1},t)1_{\\{y_{2}=0\\}}=(\omega_{0}(y_{1},0+,t)+\omega_{0}(y_{1},0-,t))1_{\\{y_{1}\geq
0,y_{2}=0\\}},$ (4.12)

as $\varepsilon\downarrow 0$ pointwise. We thus have that the integral

$\int_{D}p^{D}(0,\xi,t,y)\sigma_{\varepsilon}(\xi,0)\textrm{d}\xi\to
0\quad\textrm{ as }\varepsilon\downarrow 0,$ (4.13)

by dominated convergence as the limit of $\sigma_{\varepsilon}$ is zero almost
everywhere.

Let us consider the integral

$\int_{0}^{t}\int_{D}p^{D}(s,\xi,t,y)\rho_{\varepsilon}(\xi,s)\textrm{d}\xi\textrm{d}s.$

The following terms in (4.7)

$\phi(\pm x_{2}/\varepsilon)\left(\nu\frac{\partial^{2}\theta_{\pm}}{\partial
x_{1}^{2}}(x_{1},t)-\frac{\partial\theta_{\pm}}{\partial
t}(x_{1},t)\right)-\phi(\pm
x_{2}/\varepsilon)u^{1}(x,t)\frac{\partial\theta_{\pm}}{\partial
x_{1}}(x_{1},t),$ (4.14)

converge pointwise to

$\Big{(}\nu\frac{\partial^{2}\theta_{\pm}}{\partial
x_{1}^{2}}(x_{1},t)-\frac{\partial\theta_{\pm}}{\partial
t}(x_{1},t)\Big{)}1_{\\{x_{2}=0\\}},$ (4.15)

as $\varepsilon\downarrow 0$, due to the no-slip condition. These terms again
do not contribute to the limit as the function (4.15) is zero almost
everywhere. Similarly, the term

$\theta^{\prime}(t)\phi(|x|/\varepsilon)1_{\\{x_{1}<0\\}}$

converges pointwise to zero and does not contribute to the limit.

Consider now the following terms

$\frac{1}{\varepsilon}\phi^{\prime}(\pm
x_{2}/\varepsilon)u^{2}(x,t)\theta_{\pm}(x_{1},t),$ (4.16)

and compute the limit of these terms in weak sense. Take a smooth function
$\beta$ with compact support and write

$\displaystyle\frac{1}{\varepsilon}\int_{\mathbb{R}^{2}}\beta(x)\phi^{\prime}(x_{2}/\varepsilon)\textrm{d}x$
$\displaystyle=\int_{-\infty}^{+\infty}\int_{\frac{1}{3}}^{\frac{2}{3}}\beta(x_{1},\varepsilon
x_{2})\phi^{\prime}(x_{2})\textrm{d}x_{2}\textrm{d}x_{1}$
$\displaystyle=-\varepsilon\int_{-\infty}^{+\infty}\int_{\frac{1}{3}}^{\frac{2}{3}}\frac{\partial\beta}{\partial
x_{2}}(x_{1},\varepsilon x_{2})\phi(x_{2})\textrm{d}x_{2}\textrm{d}x_{1}$
$\displaystyle-\int_{-\infty}^{+\infty}\beta\Big{(}x_{1},\frac{\varepsilon}{3}\Big{)}\textrm{d}x_{1}.$

Therefore,

$\begin{split}\lim_{\varepsilon\to
0}\frac{1}{\varepsilon}\int_{\mathbb{R}^{2}}\beta(x)\phi^{\prime}(x_{2}/\varepsilon)\textrm{d}x=-\int_{-\infty}^{+\infty}\beta(x_{1},0)\textrm{d}x_{1},\end{split}$
(4.17)

which is zero for terms in (4.16) due to the no-slip condition. Similarly, the
integral with $\theta_{-}$ vanishes in the limit. For the terms

$\Big{(}-\frac{1}{\varepsilon}\theta(t)\phi^{\prime}(|x|/\varepsilon)\frac{u^{1}(x_{1},x_{2})x_{1}+u^{2}(x_{1},x_{2})x_{2}}{|x|}+\frac{\nu}{\varepsilon}\theta(t)\phi^{\prime}(|x|/\varepsilon)\frac{1}{|x|}\Big{)}1_{\\{x_{1}<0\\}},$
(4.18)

we write in polar coordinates

$\displaystyle\frac{1}{\varepsilon}\int_{\mathbb{R}^{2}}\beta(x)\frac{\phi^{\prime}(|x|/\varepsilon)}{|x|}\textrm{d}x$
$\displaystyle=\int_{0}^{2\pi}\int_{\frac{1}{3}}^{\frac{2}{3}}\beta(\varepsilon
r,\psi)\phi^{\prime}(r)\textrm{d}r\textrm{d}\psi$
$\displaystyle=-\varepsilon\int_{0}^{2\pi}\int_{\frac{1}{3}}^{\frac{2}{3}}\frac{\partial\beta}{\partial
r}(\varepsilon r,\psi)\phi(r)\textrm{d}r\textrm{d}\psi$
$\displaystyle-\int_{0}^{2\pi}\beta\Big{(}\frac{\varepsilon}{3},\psi\Big{)}\textrm{d}\psi.$

Hence

$\begin{split}\lim_{\varepsilon\to
0}\frac{1}{\varepsilon}\int_{\mathbb{R}^{2}}\beta(x)\frac{\phi^{\prime}(|x|/\varepsilon)}{|x|}\textrm{d}x=-\lim_{r\to
0}\int_{0}^{2\pi}\beta(r\cos\psi,r\sin\psi)\textrm{d}\psi,\end{split}$ (4.19)

which is zero for (4.18).

Finally, consider the terms of the form

$\frac{\nu}{\varepsilon^{2}}\theta_{\pm}(x_{1},t)\phi^{\prime\prime}(\pm
x_{2}/\varepsilon).$ (4.20)

Again, take a function $\beta$ and write

$\displaystyle\frac{1}{\varepsilon^{2}}\int_{\mathbb{R}^{2}}\beta(x)\phi^{\prime\prime}(x_{2}/\varepsilon)\textrm{d}x$
$\displaystyle=\frac{1}{\varepsilon}\int_{-\infty}^{+\infty}\int_{\frac{1}{3}}^{\frac{2}{3}}\beta(x_{1},\varepsilon
x_{2})\phi^{\prime\prime}(x_{2})\textrm{d}x_{2}\textrm{d}x_{1}$
$\displaystyle=-\int_{-\infty}^{+\infty}\int_{\frac{1}{3}}^{\frac{2}{3}}\frac{\partial\beta}{\partial
x_{2}}(x_{1},\varepsilon
x_{2})\phi^{\prime}(x_{2})\textrm{d}x_{2}\textrm{d}x_{1}$
$\displaystyle=\varepsilon\int_{-\infty}^{+\infty}\int_{\frac{1}{3}}^{\frac{2}{3}}\frac{\partial^{2}\beta}{\partial
x_{2}^{2}}(x_{1},\varepsilon x_{2})\phi(x_{2})\textrm{d}x_{2}\textrm{d}x_{1}$
$\displaystyle+\int_{-\infty}^{+\infty}\frac{\partial\beta}{\partial
x_{2}}\Big{(}x_{1},\frac{\varepsilon}{3}\Big{)}\textrm{d}x_{1}.$

Therefore,

$\lim_{\varepsilon\to
0}\frac{1}{\varepsilon^{2}}\int_{\mathbb{R}^{2}}\beta(x)\phi^{\prime\prime}(x_{2}/\varepsilon)\textrm{d}x=\int_{-\infty}^{+\infty}\frac{\partial\beta}{\partial
x_{2}}(x_{1},0)\textrm{d}x_{1},$ (4.21)

and the corresponding term in $\rho_{\varepsilon}$ contributes to the limit as

$\lim_{\varepsilon\to
0}\frac{\nu}{\varepsilon^{2}}\int_{D}p^{D}(s,\xi,t,y)\theta_{+}(\xi_{1},s)\phi^{\prime\prime}(\xi_{2}/\varepsilon)\textrm{d}\xi=\nu\int_{0}^{+\infty}\frac{\partial}{\partial\xi_{2}}p^{D}(s,(\xi_{1},0+),t,y)\theta_{+}(\xi_{1},s)\textrm{d}\xi_{1},$
(4.22)

and we also have a similar contribution from the integral with $\theta_{-}$.
Similarly, for the term

$\frac{\nu}{\varepsilon^{2}}\theta(t)\phi^{\prime\prime}(|x|/\varepsilon)1_{\\{x_{1}<0\\}},$
(4.23)

we write

$\displaystyle\frac{1}{\varepsilon^{2}}\int_{\mathbb{R}^{2}}\beta(x)\phi^{\prime\prime}(|x|/\varepsilon)\textrm{d}x$
$\displaystyle=\int_{0}^{2\pi}\int_{\frac{1}{3}}^{\frac{2}{3}}\beta(\varepsilon
r,\psi)\phi^{\prime\prime}(r)r\textrm{d}r\textrm{d}\psi$
$\displaystyle=-\int_{0}^{2\pi}\int_{\frac{1}{3}}^{\frac{2}{3}}\Big{(}\beta(\varepsilon
r,\psi)+\varepsilon r\frac{\partial\beta}{\partial r}(\varepsilon
r,\psi)\Big{)}\phi^{\prime}(r)\textrm{d}r\textrm{d}\psi$
$\displaystyle=\int_{0}^{2\pi}\beta\Big{(}\frac{\varepsilon}{3},\psi\Big{)}\textrm{d}\psi+\varepsilon\int_{0}^{2\pi}\int_{\frac{1}{3}}^{\frac{2}{3}}\frac{\partial\beta}{\partial
r}(\varepsilon r,\psi)\phi(r)\textrm{d}r\textrm{d}\psi$
$\displaystyle-\varepsilon\int_{0}^{2\pi}\int_{\frac{1}{3}}^{\frac{2}{3}}r\frac{\partial\beta}{\partial
r}(\varepsilon r,\psi)\phi^{\prime}(r)\textrm{d}r\textrm{d}\psi.$

Therefore,

$\begin{split}\lim_{\varepsilon\to
0}\frac{1}{\varepsilon^{2}}\int_{\mathbb{R}^{2}}\beta(x)\phi^{\prime\prime}(|x|/\varepsilon)\textrm{d}x=\lim_{r\to
0}\int_{0}^{2\pi}\beta(r\cos\psi,r\sin\psi)\textrm{d}\psi,\end{split}$ (4.24)

which is zero for the term (4.23). ∎

Therefore, the Biot-Savart law implies the following representation for the
velocity $u$.

######  Lemma 4.2.

The following integral formula holds:

$\displaystyle u^{i}(x,t)$
$\displaystyle=\int_{D}\left(\int_{D}K^{i}(x,y)p^{D}(0,\xi,t,y)\textrm{d}y\right)\omega_{0}(\xi)\textrm{d}\xi$
$\displaystyle+\int_{0}^{t}\int_{D}\left(\int_{D}K^{i}(x,y)p^{D}(s,\xi,t,y)\textrm{d}y\right)G(\xi,s)\textrm{d}\xi\textrm{d}s$
$\displaystyle+\nu\int_{0}^{t}\int_{0}^{\infty}\left(\int_{D}K^{i}(x,y)\frac{\partial}{\partial\xi_{2}}p^{D}(s,(\xi_{1},0+),t,y)\textrm{d}y\right)\theta_{+}(\xi_{1},s)\textrm{d}\xi_{1}\textrm{d}s$
$\displaystyle-\nu\int_{0}^{t}\int_{0}^{\infty}\left(\int_{D}K^{i}(x,y)\frac{\partial}{\partial\xi_{2}}p^{D}(s,(\xi_{1},0-),t,y)\textrm{d}y\right)\theta_{-}(\xi_{1},s)\textrm{d}\xi_{1}\textrm{d}s$
(4.25)

for every $x\in D$ and $t>0$.

###### Proof.

As noticed above, the formula is implied by the Biot-Savart law (3.3). Note
that the following term
$(\omega_{0}(y_{1},0+,t)+\omega_{0}(y_{1},0-,t))1_{\\{y\in\partial D\\}}$ does
not contribute to the representation as it is zero outside the boundary. ∎

Thanks to the integral representation (4.25), we may establish the random
vortex dynamics as follows.

Firstly extend $u(x,t)$ to be zero if $x\notin D$, and define the Taylor
diffusion process family $X^{\xi,s}$ as the weak solution to the stochastic
differential equation

$\begin{cases}\textrm{d}X_{t}^{\xi,s}=u(X_{t}^{\xi,s},t)\textrm{d}t+\sqrt{2\nu}\textrm{d}B_{t},&\textrm{
for }t\geq s,\\\ X_{t}^{\xi,s}=\xi,&\textrm{ for }t\leq s.\end{cases}$ (4.26)

Let

$\tau_{\xi,s}=\inf\left\\{t\geq s:X_{t}^{\xi,s}\in\partial D\right\\}.$

Then

$\int_{D}K^{i}(x,y)p^{D}(s,\xi,t,y)\textrm{d}y=\mathbb{E}\left[K^{i}(x,X_{t}^{\xi,s}):t<\tau_{\xi,s}\right],$
$\int_{D}K^{i}(x,y)\frac{\partial}{\partial\xi_{2}}p^{D}(s,(\xi_{1},0+),t,y)\textrm{d}y=\left.\frac{\partial}{\partial\xi_{2}}\right|_{\xi_{2}=0+}\mathbb{E}\left[K^{i}(x,X_{t}^{\xi,s}):t<\tau_{\xi,s}\right]$

and

$\int_{D}K^{i}(x,y)\frac{\partial}{\partial\xi_{2}}p^{D}(s,(\xi_{1},0-),t,y)\textrm{d}y=\left.\frac{\partial}{\partial\xi_{2}}\right|_{\xi_{2}=0-}\mathbb{E}\left[K^{i}(x,X_{t}^{\xi,s}):t<\tau_{\xi,s}\right].$

Therefore the velocity $u(x,t)$ can be written as the following

$\displaystyle u^{i}(x,t)$
$\displaystyle=\int_{D}\mathbb{E}\left[K^{i}(x,X_{t}^{\xi,0}):t<\tau_{\xi,0}\right]\omega_{0}(\xi)\textrm{d}\xi$
$\displaystyle+\int_{0}^{t}\int_{D}\mathbb{E}\left[K^{i}(x,X_{t}^{\xi,s}):t<\tau_{\xi,s}\right]G(\xi,s)\textrm{d}\xi\textrm{d}s$
$\displaystyle+\nu\int_{0}^{t}\int_{0}^{\infty}\left.\frac{\partial}{\partial\xi_{2}}\right|_{\xi_{2}=0+}\mathbb{E}\left[K^{i}(x,X_{t}^{\xi,s}):t<\tau_{\xi,s}\right]\theta_{+}(\xi_{1},s)\textrm{d}\xi_{1}\textrm{d}s$
$\displaystyle-\nu\int_{0}^{t}\int_{0}^{\infty}\left.\frac{\partial}{\partial\xi_{2}}\right|_{\xi_{2}=0-}\mathbb{E}\left[K^{i}(x,X_{t}^{\xi,s}):t<\tau_{\xi,s}\right]\theta_{-}(\xi_{1},s)\textrm{d}\xi_{1}\textrm{d}s.$
(4.27)

###  4.2 Flows past a wedge type obstacle

In this subsection we derive a representation similar to the above for the
wedge case. Here we will use the local coordinates for boundary components,
and this idea will be used to derive a representation for more general domains
subsequently.

Recall that the wedge obstacle is given by
$\Lambda=\left\\{x=r\textrm{e}^{\textrm{i}\theta}:r\geq 0\textrm{ and
}-\alpha\leq\theta\leq\alpha\right\\}$ for a fixed
$\alpha\in(0,\frac{\pi}{2})$, so the domain
$D=\mathbb{R}^{2}\setminus\Lambda$. Note that the boundary $\partial
D=\\{x=re^{i\theta}:r\geq 0,\theta=\pm\alpha\\}$ has two components denoted
$\partial D^{\pm}$ correspondingly. We also introduce normal coordinates
$n^{\pm},\tau^{\pm}$ in neighbourhood of $\partial D^{\pm}$ which are given by
the following transformation

$\left\\{\begin{aligned} &\tau^{\pm}=x_{1}\cos\alpha\pm x_{2}\sin\alpha,\\\
&n^{\pm}=-x_{1}\sin\alpha\pm x_{2}\cos\alpha.\end{aligned}\right.$ (4.28)

In these coordinates, we denote the components of the vector field $u$ by
$u^{\tau^{\pm}},u^{n^{\pm}}$. Note that the boundary components are given by
$\partial D^{\pm}=\\{\tau^{\pm}\geq 0,n^{\pm}=0\\}$. In this case, the no-slip
condition reads

$u(x,t)=\lim_{n^{\pm}\to 0}u(\tau^{\pm},n^{\pm},t)=0,$ (4.29)

for any $x\in\partial D^{\pm}$. As in (4.3), the stress at the boundary is
given by

$\theta_{\pm}(\tau^{\pm},t)=\begin{cases}-\left.\frac{\partial
u^{\tau^{\pm}}}{\partial
n^{\pm}}(\tau^{\pm},n^{\pm},t)\right|_{n^{\pm}=0+}&\textrm{ for
}\tau^{\pm}\geq 0,\\\ 0&\textrm{ for }\tau^{\pm}<0,\end{cases}$ (4.30)

as the tangential derivatives $\frac{\partial u}{\partial\tau^{\pm}}$ vanish
due to the no-slip condition. Note that the derivative operators are given by

$\left\\{\begin{aligned}
&\frac{\partial}{\partial\tau^{\pm}}=\cos\alpha\frac{\partial}{\partial
x_{1}}\pm\sin\alpha\frac{\partial}{\partial x_{2}},\\\
&\frac{\partial}{\partial n^{\pm}}=-\sin\alpha\frac{\partial}{\partial
x_{1}}\pm\cos\alpha\frac{\partial}{\partial x_{2}};\end{aligned}\right.$
(4.31)

due to the definition (4.28).

We introduce the extension $\sigma_{\varepsilon}$ as follows

$\sigma_{\varepsilon}(x_{1},x_{2},t)=\begin{cases}\theta_{+}(\tau^{+},t)\phi(n^{+}/\varepsilon),&\text{
if }\arg x\in[\alpha,\alpha+\frac{\pi}{2}],\\\
\theta(t)\phi(|x|/\varepsilon),&\text{ if }\arg
x\in(\alpha+\frac{\pi}{2},\frac{3\pi}{2}-\alpha),\\\
\theta_{-}(\tau^{-},t)\phi(n^{-}/\varepsilon),&\text{ if }\arg
x\in[\frac{3\pi}{2}-\alpha,2\pi-\alpha];\end{cases}$ (4.32)

where the coordinates $n^{\pm},\tau^{\pm}$ are given as functions of
$x_{1},x_{2}$ in (4.28), and $\theta(t)$ denotes the value of
$\theta_{+}(0,t)=\theta_{-}(0,t)$. Then
$W^{\varepsilon}=\omega-\sigma_{\varepsilon}$ as before satisfies the
following equation

$\left(\frac{\partial}{\partial
t}+u\cdot\nabla-\nu\Delta\right)W^{\varepsilon}=G+\rho_{\varepsilon}\quad\textrm{
in }D,$ (4.33)

with homogeneous boundary condition

$\lim_{n^{\pm}\to 0+}W^{\varepsilon}(\tau^{\pm},n^{\pm},t)=0\quad\textrm{ for
}\tau^{\pm}\geq 0,$ (4.34)

and initial data

$\displaystyle W_{0}^{\varepsilon}(x)$
$\displaystyle=\omega_{0}(x_{1},x_{2})-\theta_{+}(\tau^{+},0)\phi(n^{+}/\varepsilon)1_{\\{\arg
x\in[\alpha,\alpha+\frac{\pi}{2}]\\}}$
$\displaystyle-\theta_{-}(\tau^{-},0)\phi(n^{-}/\varepsilon)1_{\\{\arg
x\in[\frac{3\pi}{2}-\alpha,2\pi-\alpha]\\}}$
$\displaystyle-\theta(0)\phi(|x|/\varepsilon)1_{\\{\arg
x\in(\alpha+\frac{\pi}{2},\frac{3\pi}{2}-\alpha)\\}},$ (4.35)

for $x\in D$.

The function $\rho_{\varepsilon}$ is given by

$\rho_{\varepsilon}=\nu\Delta\sigma_{\varepsilon}-\frac{\partial\sigma_{\varepsilon}}{\partial
t}-u\cdot\nabla\sigma_{\varepsilon},$ (4.36)

and, using the definition of $\sigma_{\varepsilon}$, can be written as follows

$\displaystyle\rho_{\varepsilon}(x,t)$
$\displaystyle=\Big{(}\nu\frac{\partial^{2}\theta_{+}}{\partial(\tau^{+})^{2}}(\tau^{+},t)\phi(n^{+}/\varepsilon)+\frac{\nu}{\varepsilon^{2}}\theta_{+}(\tau^{+},t)\phi^{\prime\prime}(n^{+}/\varepsilon)-\frac{\partial\theta_{+}}{\partial
t}(\tau^{+},t)\phi(n^{+}/\varepsilon)$
$\displaystyle-u^{\tau^{+}}(\tau^{+},n^{+},t)\frac{\partial\theta_{+}}{\partial\tau^{+}}(\tau^{+},t)\phi(n^{+}/\varepsilon)-\frac{1}{\varepsilon}u^{n^{+}}(\tau^{+},n^{+},t)\theta_{+}(\tau^{+},t)\phi^{\prime}(n^{+}/\varepsilon)\Big{)}1_{\\{\arg
x\in[\alpha,\alpha+\frac{\pi}{2}]\\}}$
$\displaystyle+\Big{(}\nu\frac{\partial^{2}\theta_{-}}{\partial(\tau^{-})^{2}}(\tau^{-},t)\phi(n^{-}/\varepsilon)+\frac{\nu}{\varepsilon^{2}}\theta_{-}(\tau^{-},t)\phi^{\prime\prime}(n^{-}/\varepsilon)-\frac{\partial\theta_{-}}{\partial
t}(\tau^{-},t)\phi(n^{-}/\varepsilon)$
$\displaystyle-u^{\tau^{-}}(\tau^{-},n^{-},t)\frac{\partial\theta_{-}}{\partial\tau^{-}}(\tau^{-},t)\phi(n^{-}/\varepsilon)-\frac{1}{\varepsilon}u^{n^{-}}(\tau^{-},n^{-},t)\theta_{-}(\tau^{-},t)\phi^{\prime}(n^{-}/\varepsilon)\Big{)}1_{\\{\arg
x\in[\frac{3\pi}{2}-\alpha,2\pi-\alpha]\\}}$
$\displaystyle+\Big{(}-\frac{1}{\varepsilon}\theta(t)\phi^{\prime}(|x|/\varepsilon)\frac{u^{1}(x_{1},x_{2})x_{1}+u^{2}(x_{1},x_{2})x_{2}}{|x|}-\theta^{\prime}(t)\phi(|x|/\varepsilon)+\frac{\nu}{\varepsilon^{2}}\theta(t)\phi^{\prime\prime}(|x|/\varepsilon)$
$\displaystyle+\frac{\nu}{\varepsilon}\theta(t)\phi^{\prime}(|x|/\varepsilon)\frac{1}{|x|}\Big{)}1_{\\{\arg
x\in(\alpha+\frac{\pi}{2},\frac{3\pi}{2}-\alpha)\\}}.$ (4.37)

As before, we write the representation for the solution $W^{\varepsilon}$ to
(4.33) in terms of the transition probability density $p^{D}(s,\xi,t,\eta)$
which implies

$\omega(y,t)=\sigma_{\varepsilon}(y,t)+\int_{D}p^{D}(0,\xi,t,y)W^{\varepsilon}(\xi,0)\textrm{d}\xi+\int_{0}^{t}\int_{D}p^{D}(s,\xi,t,y)(G(\xi,s)+\rho_{\varepsilon}(\xi,s))\textrm{d}\xi\textrm{d}s,$
(4.38)

for $y\in D$ and $t>0$.

######  Lemma 4.3.

We have the following representation:

$\displaystyle\omega(y,t)$ $\displaystyle=\omega(y,t)1_{\\{y\in\partial
D\\}}+\int_{D}p^{D}(0,\xi,t,y)\omega_{0}(\xi)\textrm{d}\xi+\int_{0}^{t}\int_{D}p^{D}(s,\xi,t,y)G(\xi,s)\textrm{d}\xi\textrm{d}s$
$\displaystyle+\nu\int_{0}^{t}\int_{0}^{+\infty}\frac{\partial}{\partial
n^{+}}p^{D}(s,(\tau^{+},0+),t,y)\theta_{+}(\tau^{+},s)\textrm{d}\tau^{+}\textrm{d}s$
$\displaystyle+\nu\int_{0}^{t}\int_{0}^{+\infty}\frac{\partial}{\partial
n^{-}}p^{D}(s,(\tau^{-},0+),t,y)\theta_{-}(\tau^{-},s)\textrm{d}\tau^{-}\textrm{d}s,$
(4.39)

for $y\in D$ and $t>0$.

###### Proof.

To obtain the result, we again let $\varepsilon\downarrow 0$ in the
representation (4.38). First, it is easy to see that due to (4.32),

$\sigma_{\varepsilon}(y,t)\to\theta_{+}(\tau^{+},t)1_{\\{\tau^{+}\geq
0,n^{+}=0\\}}+\theta_{-}(\tau^{-},t)1_{\\{\tau^{-}\geq
0,n^{-}=0\\}}=\omega(y,t)1_{\\{y\in\partial D\\}},$ (4.40)

pointwise as $\varepsilon\to 0$. This also implies that the integral

$\int_{D}p^{D}(0,\xi,t,y)\sigma_{\varepsilon}(\xi,0)\textrm{d}\xi\to
0,\quad\textrm{ as }\varepsilon\downarrow 0,$ (4.41)

by dominated convergence.

Now let us take the limit in the last integral in (4.38). It is easy to see
that in (4.2), the following terms

$\phi(n^{\pm}/\varepsilon)\left(\nu\frac{\partial^{2}\theta_{\pm}}{\partial(\tau^{\pm})^{2}}(\tau^{\pm},t)-\frac{\partial\theta_{\pm}}{\partial
t}(\tau^{\pm},t)\right)-u^{\tau^{\pm}}(\tau^{\pm},n^{\pm},t)\frac{\partial\theta_{\pm}}{\partial\tau^{\pm}}(\tau^{\pm},t)\phi(n^{\pm}/\varepsilon),$
(4.42)

as $\varepsilon\to 0$, converge pointwise to

$\Big{(}\nu\frac{\partial^{2}\theta_{\pm}}{\partial(\tau^{\pm})^{2}}(\tau^{\pm},t)-\frac{\partial\theta_{\pm}}{\partial
t}(\tau^{\pm},t)\Big{)}1_{\\{\tau^{\pm}\geq 0,n^{\pm}=0\\}},$ (4.43)

since $u|_{\partial D^{\pm}}=0$ and therefore, these terms do not contribute
to the limit. Similarly,

$-\theta^{\prime}(t)\phi(|x|/\varepsilon)1_{\\{\arg
x\in(\alpha+\frac{\pi}{2},\frac{3\pi}{2}-\alpha)\\}}$

converges to zero pointwise and does not contribute to the limiting
representation.

Let us consider now the term

$\frac{1}{\varepsilon}u^{n^{+}}(\tau^{+},n^{+},t)\theta_{+}(\tau^{+},t)\phi^{\prime}(n^{+}/\varepsilon)1_{\\{\arg
x\in[\alpha,\alpha+\frac{\pi}{2}]\\}}.$ (4.44)

Take a smooth function $\beta$ with compact support and consider first the
integral

$\displaystyle\frac{1}{\varepsilon}\int_{\mathbb{R}^{2}}\beta(x)\phi^{\prime}(n^{+}/\varepsilon)1_{\\{\arg
x\in[\alpha,\alpha+\frac{\pi}{2}]\\}}\textrm{d}x$
$\displaystyle=\frac{1}{\varepsilon}\int_{0}^{+\infty}\int_{0}^{+\infty}\beta(\tau^{+},n^{+})\phi^{\prime}(n^{+}/\varepsilon)\textrm{d}n^{+}\textrm{d}\tau^{+}$
$\displaystyle=\int_{0}^{+\infty}\int_{\frac{1}{3}}^{\frac{2}{3}}\beta(\tau^{+},\varepsilon
n^{+})\phi^{\prime}(n^{+})\textrm{d}n^{+}\textrm{d}\tau^{+}$
$\displaystyle=-\int_{0}^{+\infty}\beta\Big{(}\tau^{+},\frac{\varepsilon}{3}\Big{)}\textrm{d}\tau^{+}$
$\displaystyle-\varepsilon\int_{0}^{+\infty}\int_{\frac{1}{3}}^{\frac{2}{3}}\frac{\partial\beta}{\partial
n^{+}}(\tau^{+},\varepsilon
n^{+})\phi(n^{+})\textrm{d}n^{+}\textrm{d}\tau^{+}.$

Therefore, we have

$\begin{split}\lim_{\varepsilon\to
0}\frac{1}{\varepsilon}\int_{\mathbb{R}^{2}}\beta(x)\phi^{\prime}(n^{+}/\varepsilon)1_{\\{\arg
x\in[\alpha,\alpha+\frac{\pi}{2}]\\}}\textrm{d}x=-\int_{0}^{+\infty}\beta(\tau^{+},0)\textrm{d}\tau^{+},\end{split}$
(4.45)

which is zero for terms in (4.44) due to the no-slip condition. Similarly, the
integral with the term

$\frac{1}{\varepsilon}u^{n^{-}}(\tau^{-},n^{-},t)\theta_{-}(\tau^{-},t)\phi^{\prime}(n^{-}/\varepsilon)1_{\\{\arg
x\in[\frac{3\pi}{2}-\alpha,2\pi-\alpha]\\}}.$

vanishes in the limit. Consider now the following terms

$\Big{(}-\frac{1}{\varepsilon}\theta(t)\phi^{\prime}(|x|/\varepsilon)\frac{u^{1}(x_{1},x_{2})x_{1}+u^{2}(x_{1},x_{2})x_{2}}{|x|}+\frac{\nu}{\varepsilon}\theta(t)\phi^{\prime}(|x|/\varepsilon)\frac{1}{|x|}\Big{)}1_{\\{\arg
x\in(\alpha+\frac{\pi}{2},\frac{3\pi}{2}-\alpha)\\}}.$ (4.46)

Recall that in the proof of Lemma 4.1 we had

$\begin{split}\lim_{\varepsilon\to
0}\frac{1}{\varepsilon}\int_{\mathbb{R}^{2}}\beta(x)\frac{\phi^{\prime}(|x|/\varepsilon)}{|x|}\textrm{d}x=-\lim_{r\to
0}\int_{0}^{2\pi}\beta(r\cos\psi,r\sin\psi)\textrm{d}\psi,\end{split}$ (4.47)

which is zero for the terms (4.46).

Finally, consider the terms of the form

$\frac{\nu}{\varepsilon^{2}}\theta_{\pm}(\tau^{\pm},t)\phi^{\prime\prime}(n^{\pm}/\varepsilon).$
(4.48)

Again, take a function $\beta$ and write

$\displaystyle\frac{1}{\varepsilon^{2}}\int_{\mathbb{R}^{2}}\beta(x)\phi^{\prime\prime}(n^{+}/\varepsilon)1_{\\{\arg
x\in[\alpha,\alpha+\frac{\pi}{2}]\\}}\textrm{d}x$
$\displaystyle=\frac{1}{\varepsilon^{2}}\int_{0}^{+\infty}\int_{0}^{+\infty}\beta(\tau^{+},n^{+})\phi^{\prime\prime}(n^{+}/\varepsilon)\textrm{d}n^{+}\textrm{d}\tau^{+}$
$\displaystyle=\frac{1}{\varepsilon}\int_{0}^{+\infty}\int_{\frac{1}{3}}^{\frac{2}{3}}\beta(\tau^{+},\varepsilon
n^{+})\phi^{\prime\prime}(n^{+})\textrm{d}n^{+}\textrm{d}\tau^{+}$
$\displaystyle=-\int_{0}^{+\infty}\int_{\frac{1}{3}}^{\frac{2}{3}}\frac{\partial\beta}{\partial
n^{+}}(\tau^{+},\varepsilon
n^{+})\phi^{\prime}(n^{+})\textrm{d}n^{+}\textrm{d}\tau^{+}$
$\displaystyle=\int_{0}^{+\infty}\frac{\partial\beta}{\partial
n^{+}}\Big{(}\tau^{+},\frac{\varepsilon}{3}\Big{)}\textrm{d}\tau^{+}$
$\displaystyle+\varepsilon\int_{0}^{+\infty}\int_{\frac{1}{3}}^{\frac{2}{3}}\frac{\partial^{2}\beta}{\partial(n^{+})^{2}}(\tau^{+},\varepsilon
n^{+})\phi(n^{+})\textrm{d}n^{+}\textrm{d}\tau^{+}.$

Therefore,

$\lim_{\varepsilon\to
0}\frac{1}{\varepsilon^{2}}\int_{\mathbb{R}^{2}}\beta(x)\phi^{\prime\prime}(n^{+}/\varepsilon)1_{\\{\arg
x\in[\alpha,\alpha+\frac{\pi}{2}]\\}}\textrm{d}x=\int_{0}^{+\infty}\frac{\partial\beta}{\partial
n^{+}}(\tau^{+},0)\textrm{d}\tau^{+},$ (4.49)

and the corresponding term in $\rho_{\varepsilon}$ contributes to the limit as

$\begin{split}\lim_{\varepsilon\to
0}\frac{\nu}{\varepsilon^{2}}\int_{D}p^{D}(s,\xi,t,y)\theta_{+}(\tau^{+},s)\phi^{\prime\prime}(n^{+}/\varepsilon)1_{\\{y_{2}\geq
0\\}}\textrm{d}\xi\\\ =\nu\int_{0}^{+\infty}\frac{\partial}{\partial
n^{+}}p^{D}(s,(\tau^{+},0+),t,y)\theta_{+}(\tau^{+},s)\textrm{d}\tau^{+},\end{split}$
(4.50)

and we also have a similar contribution from the integral with $\theta_{-}$.
Lastly, consider the term

$\frac{\nu}{\varepsilon^{2}}\theta(t)\phi^{\prime\prime}(|x|/\varepsilon)1_{\\{\arg
x\in(\alpha+\frac{\pi}{2},\frac{3\pi}{2}-\alpha)\\}}.$ (4.51)

Recall that in Lemma 4.1, it was shown that

$\begin{split}\lim_{\varepsilon\to
0}\frac{1}{\varepsilon^{2}}\int_{\mathbb{R}^{2}}\beta(x)\phi^{\prime\prime}(|x|/\varepsilon)\textrm{d}x=\lim_{r\to
0}\int_{0}^{2\pi}\beta(r\cos\psi,r\sin\psi)\textrm{d}\psi,\end{split}$ (4.52)

which is zero for the term (4.51). ∎

We again have the following representation for the velocity due to the Biot-
Savart law.

######  Lemma 4.4.

The following representation holds:

$\displaystyle u^{i}(x,t)$
$\displaystyle=\int_{D}\left(\int_{D}K^{i}(x,y)p^{D}(0,\xi,t,y)\textrm{d}y\right)\omega_{0}(\xi)\textrm{d}\xi$
$\displaystyle+\int_{0}^{t}\int_{D}\left(\int_{D}K^{i}(x,y)p^{D}(s,\xi,t,y)\textrm{d}y\right)G(\xi,s)\textrm{d}\xi\textrm{d}s$
$\displaystyle+\nu\int_{0}^{t}\int_{0}^{\infty}\left(\int_{D}K^{i}(x,y)\frac{\partial}{\partial
n^{+}}p^{D}(s,(\tau^{+},0+),t,y)\textrm{d}y\right)\theta_{+}(\tau^{+},s)\textrm{d}\tau^{+}\textrm{d}s$
$\displaystyle+\nu\int_{0}^{t}\int_{0}^{\infty}\left(\int_{D}K^{i}(x,y)\frac{\partial}{\partial
n^{-}}p^{D}(s,(\tau^{-},0+),t,y)\textrm{d}y\right)\theta_{-}(\tau^{-},s)\textrm{d}\tau^{-}\textrm{d}s,$
(4.53)

for any $x\in D$ and $t>0$.

We can use the Taylor diffusion processes $X^{\xi,s}$ defined in (4.26) to
represent $u$ as follows. Notice that

$\int_{D}K^{i}(x,y)p^{D}(s,\xi,t,y)\textrm{d}y=\mathbb{E}\left[K^{i}(x,X_{t}^{\xi,s}):t<\tau_{\xi,s}\right],$

and

$\int_{D}K^{i}(x,y)\left.\frac{\partial}{\partial
n^{\pm}}\right|_{n^{\pm}=0+}p^{D}(s,(\tau^{\pm},n^{\pm}),t,y)\textrm{d}y=\left.\frac{\partial}{\partial
n^{\pm}}\right|_{n^{\pm}=0+}\mathbb{E}\left[K^{i}(x,X_{t}^{\xi,s}):t<\tau_{\xi,s}\right],$

where $\tau_{\xi,s}$ is again the boundary hitting time for $X^{\xi,s}$. This
implies that

$\displaystyle
u^{i}(x,t)=\int_{D}\mathbb{E}\left[K^{i}(x,X_{t}^{\xi,0}):t<\tau_{\xi,0}\right]\omega_{0}(\xi)\textrm{d}\xi+\int_{0}^{t}\int_{D}\mathbb{E}\left[K^{i}(x,X_{t}^{\xi,s}):t<\tau_{\xi,s}\right]g(\xi,s)\textrm{d}\xi\textrm{d}s$
$\displaystyle+\nu\int_{0}^{t}\int_{0}^{+\infty}\frac{\partial}{\partial
n^{+}}\Big{|}_{n^{+}=0+}\mathbb{E}\left[K^{i}(x,X_{t}^{\xi,s}):t<\tau_{\xi,s}\right]\theta_{+}(\tau^{+},s)\textrm{d}\tau^{+}\textrm{d}s$
$\displaystyle+\nu\int_{0}^{t}\int_{0}^{+\infty}\frac{\partial}{\partial
n^{-}}\Big{|}_{n^{-}=0+}\mathbb{E}\left[K^{i}(x,X_{t}^{\xi,s}):t<\tau_{\xi,s}\right]\theta_{-}(\tau^{-},s)\textrm{d}\tau^{-}\textrm{d}s.$
(4.54)

###  4.3 General case

We consider a flow in $D\subset\mathbb{R}^{2}$ assuming it is a proper simply
connected domain. Let us fix a conformal mapping $T:D\rightarrow\mathbb{H}$.
Notice that in this case we can parameterise the domain $D$ using the
coordinates $(z_{1},z_{2})\in\mathbb{H}$, i.e. the transform is given by
$x=T^{-1}(z)$ for $x\in\mathbb{H}$. The boundary in this coordinates is given
by $\partial D=\\{z_{2}=0\\}$, and the no-slip condition is written as

$\lim_{z_{2}\to 0+}u(z_{1},z_{2},t)=0.$ (4.55)

Notice that as the velocity $u$ is a vector field, its components
$u^{z_{1}},u^{z_{2}}$ in new coordinates are transformed correspondingly. We
again introduce the boundary stress as

$\theta(z_{1},t)=-\left.\frac{\partial u^{z_{1}}}{\partial
z_{2}}(z_{1},z_{2},t)\right|_{z_{2}=0+}.$ (4.56)

Note that in the subsequent argument we have to assume that the boundary
$\partial D$ is sufficiently regular, in particular, $T$ gives a smooth
parametrisation of the boundary (except, maybe, one point in $\partial D$
which is sent to infinity by $T$). However, as for the plate and wedge cases
we were able to take care of singularities at the boundary $\partial D$, we
expect a similar formula to hold more generally, say, for domains with
piecewise smooth boundary.

Let us define the extension $\sigma_{\varepsilon}$ in new coordinates by

$\sigma_{\varepsilon}(z_{1},z_{2},t)=\theta(z_{1},t)\phi(z_{2}/\varepsilon),$
(4.57)

and introduce $W^{\varepsilon}=\omega-\sigma_{\varepsilon}$. It has to satisfy

$\left(\frac{\partial}{\partial
t}+u\cdot\nabla-\nu\Delta\right)W^{\varepsilon}=G+\rho_{\varepsilon}\quad\textrm{
in }D,$ (4.58)

with homogeneous boundary condition $\left.W^{\varepsilon}\right|_{\partial
D}=0$ written in new coordinates as

$\lim_{z_{2}\to 0+}W^{\varepsilon}(z_{1},z_{2},t)=0$

for all $z_{1}$. The initial data for $W^{\varepsilon}$ is

$W_{0}^{\varepsilon}(z)=\omega_{0}(z_{1},z_{2})-\theta(z_{1},0)\phi(z_{2}/\varepsilon),$
(4.59)

and $\rho_{\varepsilon}$ is given by

$\rho_{\varepsilon}=\nu\Delta\sigma_{\varepsilon}-\frac{\partial\sigma_{\varepsilon}}{\partial
t}-u\cdot\nabla\sigma_{\varepsilon}.$ (4.60)

Therefore, in new coordinates

$\displaystyle\rho_{\varepsilon}(z,t)=-\frac{\partial\theta}{\partial
t}(z_{1},t)\phi(z_{2}/\varepsilon)$
$\displaystyle-u^{z_{1}}(z_{1},z_{2},t)\frac{\partial\theta}{\partial
z_{1}}(z_{1},t)\phi(z_{2}/\varepsilon)-\frac{1}{\varepsilon}u^{z_{2}}(z_{1},z_{2},t)\theta(z_{1},t)\phi^{\prime}(z_{2}/\varepsilon)$
$\displaystyle+\nu|T^{\prime}|^{2}\Big{(}\frac{\partial^{2}\theta}{\partial
z_{1}^{2}}(z_{1},t)\phi(z_{2}/\varepsilon)+\frac{1}{\varepsilon^{2}}\theta(z_{1},t)\phi^{\prime\prime}(z_{2}/\varepsilon)\Big{)},$
(4.61)

where $T^{\prime}$ denotes the derivative of $T$ as a complex function.

Writing the solution to (4.58) in terms of the transition density
$p^{D}(s,\xi,t,\eta)$ we have

$\omega(y,t)=\sigma_{\varepsilon}(y,t)+\int_{D}p^{D}(0,\xi,t,y)W^{\varepsilon}_{0}(\xi)\textrm{d}\xi+\int_{0}^{t}\int_{D}p^{D}(s,\xi,t,y)(G(\xi,s)+\rho_{\varepsilon}(\xi,s))\textrm{d}\xi\textrm{d}s,$
(4.62)

for $y\in D$ and $t>0$.

######  Lemma 4.5.

We have the following representation:

$\displaystyle\omega(y,t)=\omega(y,t)1_{\\{y\in\partial D\\}}$
$\displaystyle+\int_{D}p^{D}(0,\xi,t,y)\omega_{0}(\xi)\textrm{d}\xi+\int_{0}^{t}\int_{D}p^{D}(s,\xi,t,y)G(\xi,s)\textrm{d}\xi\textrm{d}s$
$\displaystyle+\nu\int_{0}^{t}\int_{-\infty}^{+\infty}\frac{\partial}{\partial
z_{2}}\Big{|}_{z_{2}=0}p^{D}(s,z,t,y)\theta(z_{1},s)\textrm{d}z_{1}\textrm{d}s,$
(4.63)

for $y\in D$ and $t>0$.

###### Proof.

We take as usual the limit in (4.62) as $\varepsilon\to 0$. Notice that as
$\varepsilon\to 0$,

$\sigma_{\varepsilon}(z,t)\to\theta(z_{1},t)1_{\\{z_{2}=0\\}},$ (4.64)

pointwise, which can be written as $\omega(y,t)1_{\\{y\in\partial D\\}}$. This
also implies that the integral

$\int_{D}p^{D}(0,\xi,t,y)\sigma_{\varepsilon}(\xi,0)\textrm{d}\xi\to
0,\quad\textrm{ as }\varepsilon\downarrow 0,$ (4.65)

by dominated convergence.

Now let us find the limit of

$\int_{0}^{t}\int_{D}p^{D}(s,\xi,t,y)\rho_{\varepsilon}(\xi,s)\textrm{d}\xi\textrm{d}s,$
(4.66)

as $\varepsilon\to 0$. Notice that the terms

$-\frac{\partial\theta}{\partial
t}(z_{1},t)\phi(z_{2}/\varepsilon)-u^{z_{1}}(z_{1},z_{2},t)\frac{\partial\theta}{\partial
z_{1}}(z_{1},t)\phi(z_{2}/\varepsilon)+\nu|T^{\prime}|^{2}\frac{\partial^{2}\theta}{\partial
z_{1}^{2}}(z_{1},t)\phi(z_{2}/\varepsilon),$ (4.67)

converge pointwise to

$\Big{(}-\frac{\partial\theta}{\partial
t}(z_{1},t)+\nu|T^{\prime}|^{2}\frac{\partial^{2}\theta}{\partial
z_{1}^{2}}(z_{1},t)\Big{)}1_{\\{z_{2}=0\\}},$ (4.68)

as $\left.u\right|_{\partial D}=0$. These terms do not contribute to the limit
as (4.68) is zero almost everywhere.

Consider now the term

$-\frac{1}{\varepsilon}u^{z_{2}}(z_{1},z_{2},t)\theta(z_{1},t)\phi^{\prime}(z_{2}/\varepsilon).$
(4.69)

Take a smooth function $\beta$ with compact support and write

$\frac{1}{\varepsilon}\int_{\mathbb{R}^{2}}\beta(x)\phi^{\prime}(z_{2}(x)/\varepsilon)\textrm{d}x=\frac{1}{\varepsilon}\int_{-\infty}^{+\infty}\int_{0}^{+\infty}\beta(z_{1},z_{2})\phi^{\prime}(z_{2}/\varepsilon)J(T^{-1})\textrm{d}z_{2}\textrm{d}z_{1},$

where $J(T^{-1})$ is the Jacobian of the transform $x=T^{-1}(z)$. Write the
latter integral as

$\displaystyle\int_{-\infty}^{+\infty}\int_{\frac{1}{3}}^{\frac{2}{3}}\beta(z_{1},\varepsilon
z_{2})\phi^{\prime}(z_{2})J(T^{-1})(z_{1},\varepsilon
z_{2})\textrm{d}z_{2}\textrm{d}z_{1}$
$\displaystyle=-\int_{-\infty}^{+\infty}\beta\left(z_{1},\frac{\varepsilon}{3}\right)J(T^{-1})\left(z_{1},\frac{\varepsilon}{3}\right)\textrm{d}z_{1}$
$\displaystyle-\varepsilon\int_{-\infty}^{+\infty}\int_{\frac{1}{3}}^{\frac{2}{3}}\phi(z_{2})\frac{\partial\beta}{\partial
z_{2}}(z_{1},\varepsilon z_{2})J(T^{-1})(z_{1},\varepsilon
z_{2})\textrm{d}z_{2}\textrm{d}z_{1}$
$\displaystyle-\varepsilon\int_{-\infty}^{+\infty}\int_{\frac{1}{3}}^{\frac{2}{3}}\phi(z_{2})\beta(z_{1},\varepsilon
z_{2})\frac{\partial}{\partial z_{2}}J(T^{-1})(z_{1},\varepsilon
z_{2})\textrm{d}z_{2}\textrm{d}z_{1}.$

Therefore, the limit of the term (4.69) is equal to zero due to the no-slip
condition as

$\begin{split}\lim_{\varepsilon\to
0}\frac{1}{\varepsilon}\int_{\mathbb{R}^{2}}\beta(x)\phi^{\prime}(z_{2}(x)/\varepsilon)\textrm{d}x=-\int_{-\infty}^{+\infty}\beta(z_{1},0)J(T^{-1})(z_{1},0)\textrm{d}z_{1}.\end{split}$
(4.70)

Let us consider the last term

$\frac{\nu}{\varepsilon^{2}}|T^{\prime}|^{2}\theta(z_{1},t)\phi^{\prime\prime}(z_{2}/\varepsilon).$
(4.71)

Again, taking a smooth compactly supported function $\beta$, write

$\frac{1}{\varepsilon^{2}}\int_{\mathbb{R}^{2}}\beta(x)|T^{\prime}|^{2}\phi^{\prime\prime}(z_{2}(x)/\varepsilon)\textrm{d}x=\frac{1}{\varepsilon^{2}}\int_{-\infty}^{+\infty}\int_{0}^{+\infty}\beta(z_{1},z_{2})\phi^{\prime\prime}(z_{2}/\varepsilon)\textrm{d}z_{2}\textrm{d}z_{1},$

since the Jacobian $J(T)=|T^{\prime}|^{2}$. Therefore,

$\displaystyle\frac{1}{\varepsilon}\int_{-\infty}^{+\infty}\int_{\frac{1}{3}}^{\frac{2}{3}}\beta(z_{1},\varepsilon
z_{2})\phi^{\prime\prime}(z_{2})\textrm{d}z_{2}\textrm{d}z_{1}$
$\displaystyle=-\int_{-\infty}^{+\infty}\int_{\frac{1}{3}}^{\frac{2}{3}}\frac{\partial\beta}{\partial
z_{2}}(z_{1},\varepsilon
z_{2})\phi^{\prime}(z_{2})\textrm{d}z_{2}\textrm{d}z_{1}$
$\displaystyle=\int_{-\infty}^{+\infty}\frac{\partial\beta}{\partial
z_{2}}\Big{(}z_{1},\frac{\varepsilon}{3}\Big{)}\textrm{d}z_{1}$
$\displaystyle+\varepsilon\int_{-\infty}^{+\infty}\int_{\frac{1}{3}}^{\frac{2}{3}}\frac{\partial^{2}\beta}{\partial
z_{2}^{2}}(z_{1},\varepsilon z_{2})\phi(z_{2})\textrm{d}z_{2}\textrm{d}z_{1},$

and the limit

$\lim_{\varepsilon\to
0}\frac{1}{\varepsilon^{2}}\int_{\mathbb{R}^{2}}\beta(x)|T^{\prime}|^{2}\phi^{\prime\prime}(z_{2}(x)/\varepsilon)\textrm{d}x=\int_{-\infty}^{+\infty}\frac{\partial\beta}{\partial
z_{2}}(z_{1},0)\textrm{d}z_{1}.$ (4.72)

Then it is the only term contributing to the limit, and we have that

$\int_{0}^{t}\int_{D}p^{D}(s,\xi,t,y)\rho_{\varepsilon}(\xi,s)\textrm{d}\xi\textrm{d}s\to\nu\int_{0}^{t}\int_{-\infty}^{+\infty}\frac{\partial}{\partial
z_{2}}\Big{|}_{z_{2}=0+}p^{D}(s,z,t,y)\theta(z_{1},s)\textrm{d}z_{1}\textrm{d}s,$

as $\varepsilon\to 0$. ∎

######  Lemma 4.6.

We have the following:

$\displaystyle u^{i}(x,t)$
$\displaystyle=\int_{D}\left(\int_{D}K^{i}(x,y)p^{D}(0,\xi,t,y)\textrm{d}y\right)\omega_{0}(\xi)\textrm{d}\xi$
$\displaystyle+\int_{0}^{t}\int_{D}\left(\int_{D}K^{i}(x,y)p^{D}(s,\xi,t,y)\textrm{d}y\right)G(\xi,s)\textrm{d}\xi\textrm{d}s$
$\displaystyle+\nu\int_{0}^{t}\int_{-\infty}^{+\infty}\left(\int_{D}K^{i}(x,y)\frac{\partial}{\partial
z_{2}}\Big{|}_{z_{2}=0+}p^{D}(s,(z_{1},z_{2}),t,y)\textrm{d}y\right)\theta(z_{1},s)\textrm{d}z_{1}\textrm{d}s,$
(4.73)

for any $x\in D$ and $t>0$ — note that the last integral is written in
coordinates $z_{1},z_{2}$.

The lemma above follows as usual from the Biot-Savart law. Moreover, we can
again write the representation for the velocity $u$ in terms of the Taylor
diffusions $X^{\xi,s}$ given in (4.26) as follows

$\displaystyle
u^{i}(x,t)=\int_{D}\mathbb{E}\left[K^{i}(x,X_{t}^{\xi,0}):t<\tau_{\xi,0}\right]\omega_{0}(\xi)\textrm{d}\xi+\int_{0}^{t}\int_{D}\mathbb{E}\left[K^{i}(x,X_{t}^{\xi,s}):t<\tau_{\xi,s}\right]G(\xi,s)\textrm{d}\xi\textrm{d}s$
$\displaystyle+\nu\int_{0}^{t}\int_{-\infty}^{+\infty}\frac{\partial}{\partial
z_{2}}\Big{|}_{z_{2}=0+}\mathbb{E}\left[K^{i}(x,X_{t}^{(z_{1},z_{2}),s}):t<\tau_{z,s}\right]\theta(z_{1},s)\textrm{d}z_{1}\textrm{d}s.$
(4.74)

##  5 Monte-Carlo simulations

###  5.1 Numerical scheme

The integral representations for the velocity $u$ we have obtained so far,
together with the Taylor diffusions, form a closed system. Therefore, we can
discretise the corresponding representations to establish the following
numerical schemes. For our simulations we use the representation for $u$ as in
Theorem 2.4 since in this case we work only with the diffusion processes
$X_{t}^{\eta}$ initialised at time $t=0$ (compared to (4.1), (4.2)). Thus we
use the representation

$\displaystyle u^{i}(x,t)$
$\displaystyle=\int_{D}K^{i}(x,\eta)\sigma_{\varepsilon}(\eta,t)\textrm{d}\eta+\int_{D}\mathbb{E}\left[K^{i}(x,X_{t}^{\eta})1_{\left\\{t<\zeta(X^{\eta}\circ\tau_{t})\right\\}}\right]W_{\varepsilon}(\eta,0)\textrm{d}\eta$
$\displaystyle+\int_{0}^{t}\int_{D}\mathbb{E}\left[1_{\\{t-s<\zeta(X^{\eta}\circ\tau_{t})\\}}K^{i}(x,X_{t}^{\eta})g_{\varepsilon}(X_{s}^{\eta},s)\right]\textrm{d}\eta\textrm{d}s,$
(5.1)

for the domains $D$ in Subsections 4.1 and 4.2. In fact the derivations in
these sections imply that we can pass to the limit as $\varepsilon\to 0$ in
every term except for the terms with
$\theta_{\pm}(\cdot\,,t)\phi^{\prime\prime}(\cdot\,/\varepsilon)$ for which we
might not have the same representation in the limit. We therefore approximate
the above representation by the following

$\displaystyle u^{i}(x,t)$
$\displaystyle=\int_{D}\mathbb{E}\left[1_{\left\\{t<\zeta(X^{\eta}\circ\tau_{t})\right\\}}K^{i}(x,X_{t}^{\xi})\right]\omega_{0}(\xi)\textrm{d}\xi+\int_{0}^{t}\int_{D}\mathbb{E}\left[1_{\\{t-s<\zeta(X^{\xi}\circ\tau_{t})\\}}K^{i}(x,X_{t}^{\xi})G(X_{s}^{\xi},s)\right]\textrm{d}\xi\textrm{d}s$
$\displaystyle+\frac{\nu}{\varepsilon^{2}}\int_{0}^{t}\int_{D}\mathbb{E}\left[1_{\\{t-s<\zeta(X^{\xi}\circ\tau_{t})\\}}K^{i}(x,X_{t}^{\xi})\theta_{+}(X_{s}^{\xi},s)\phi^{\prime\prime}(X_{s}^{\xi}/\varepsilon)\right]\textrm{d}\xi\textrm{d}s$
$\displaystyle+\frac{\nu}{\varepsilon^{2}}\int_{0}^{t}\int_{D}\mathbb{E}\left[1_{\\{t-s<\zeta(X^{\xi}\circ\tau_{t})\\}}K^{i}(x,X_{t}^{\xi})\theta_{-}(X_{s}^{\xi},s)\phi^{\prime\prime}(X_{s}^{\xi}/\varepsilon)\right]\textrm{d}\xi\textrm{d}s,$
(5.2)

with some small $\varepsilon$. Notice that for our choice of the cutoff
function $\phi$ as in (4.11), the second derivative $\phi^{\prime\prime}$ is
supported on $[1/3,2/3]$ and therefore the last two integrals are taken over a
layer close to the boundary.

Since the kernel $K(x,y)$ is singular, we have to mollify it, e.g. using the
cutoff function $f_{\delta}(y)=(1-\exp(-|y|/\delta))$. Notice that in the
general expression for $K$ given in Lemma 3.1 we also have the derivatives
$\dfrac{\partial T^{i}}{\partial y_{j}}$ which are singular in our case
according to Lemmas 3.3 and 3.4. We therefore denote $K_{\delta}(x,y)$ the
kernel obtained after desingularisation of $k^{\pm}(x,y)$ and $\dfrac{\partial
T^{i}}{\partial y_{j}}$.

We first have to discretise the domain $D$ choosing lattices to represent
points in $D$.

Flat-plate. For the domain $D=\mathbb{R}^{2}\setminus\\{(x_{1},0):x_{1}\geq
0\\}$, we introduce the following lattices.

  1. 1.

For each boundary component $\partial D^{\pm}$, we define the corresponding
thin boundary layer lattice $D_{b}^{\pm}$. Choose mesh sizes $h_{1},h_{2}$ and
numbers of points $N_{1},N_{2}$, and define the following lattice points

$x_{b\pm}^{i_{1}i_{2}}=(i_{1}h_{1},\pm i_{2}h_{2}),$ (5.3)

for $0\leq i_{1}\leq N_{1}$ and $0\leq i_{2}\leq N_{2}$. The number of points
used in these layers is $2(N_{1}+1)(N_{2}+1)$.

  2. 2.

Outside the thin boundary layers, we introduce the outer layer lattice
$D_{o}$. Choose a mesh size $h_{0}$ and a number $N_{0}$, and define

$x_{o}^{i_{1}i_{2}}=(i_{1}h_{0},i_{2}h_{0}),$ (5.4)

where $-N_{0}\leq i_{1}<0,|i_{2}|\leq N_{0}$ and $0\leq i_{1}\leq
N_{0},0<|i_{2}|\leq N_{0}$. The number of points used in this layer is
$\leq(2N_{0}+1)^{2}$.

Wedge. Let us also introduce the lattice points for the domain $D=\\{x:\arg
x\in(\alpha,2\pi-\alpha)\\}$.

  1. 1.

To introduce two thin boundary layer lattices $D_{b}^{\pm}$ for the
corresponding boundary components $\partial D^{\pm}$, we choose mesh sizes
$h_{1},h_{2}$ and numbers of points $N_{1},N_{2}$. Define the following
lattice points

$x_{b\pm}^{i_{1}i_{2}}=(i_{1}h_{1}\cos\alpha-i_{2}h_{2}\sin\alpha,\pm
i_{1}h_{1}\sin\alpha\pm i_{2}h_{2}\cos\alpha),$ (5.5)

for $0\leq i_{1}\leq N_{1}$ and $0\leq i_{2}\leq N_{2}$. The number of points
used in these layers is $2(N_{1}+1)(N_{2}+1)$.

  2. 2.

We also introduce the outer layer lattice $D_{o}$ choosing a mesh size $h_{0}$
and a number $N_{0}$. We define the points

$x_{o}^{i_{1}i_{2}}=(i_{1}h_{0},i_{2}h_{0}),$ (5.6)

for $|i_{1}|,|i_{2}|\leq N_{0}$. This is a square lattice in which the domain
$D$ is given by the condition: $x_{o}^{i_{1}i_{2}}\in D$ if and only if
$i_{1}\leq 0$ or $\arg(i_{1},i_{2})\in(\alpha,2\pi-\alpha)$. The number of
points used in this layer to represent points in $D$ is $\leq(2N_{0}+1)^{2}$.

We choose a time mesh size $h$ and denote $t_{k}=kh$ for $k\geq 0$. We
initialise the processes $X^{i_{1},i_{2}}_{b\pm;t_{0}}=x_{b\pm}^{i_{1}i_{2}}$
and $X^{i_{1},i_{2}}_{o;t_{0}}=x_{o}^{i_{1}i_{2}}$ and update them for $k\geq
0$ according to

$X^{i_{1},i_{2}}_{t_{k+1}}=X^{i_{1},i_{2}}_{t_{k}}+h\hat{u}(X^{i_{1},i_{2}}_{t_{k}},t_{k})+\sqrt{2\nu}(B_{t_{k+1}}-B_{t_{k}}),$
(5.7)

where we suppressed the subscripts $o$ and $b\pm$ for uniform notation. The
drift $\hat{u}$ is given by

$\displaystyle\hat{u}(x,t_{k+1})$ $\displaystyle=\sum_{(i_{1},i_{2})\in
D}A_{i_{1},i_{2}}\omega_{i_{1},i_{2}}\mathbb{E}\left[1_{\left\\{t_{k}<\zeta(X^{i_{1},i_{2}}\circ\tau_{t_{k}})\right\\}}K_{\delta}(x,X^{i_{1},i_{2}}_{t_{k}})\right]$
$\displaystyle+\sum_{(i_{1},i_{2})\in
D}A_{i_{1},i_{2}}\sum_{l=0}^{k}hG_{i_{1},i_{2};t_{l}}\mathbb{E}\left[1_{\\{t_{k}-t_{l}<\zeta(X^{i_{1},i_{2}}\circ\tau_{t_{k}})\\}}K_{\delta}(x,X^{i_{1},i_{2}}_{t_{l}})\right]$
$\displaystyle+\frac{\nu}{\varepsilon^{2}}\sum_{(i_{1},i_{2})\in
D_{b}^{-}}h_{1}h_{2}\sum_{l=0}^{k}h\mathbb{E}\left[1_{\\{t_{k}-t_{l}<\zeta(X^{i_{1},i_{2}}\circ\tau_{t_{k}})\\}}K_{\delta}(x,X^{i_{1},i_{2}}_{t_{k}})\theta_{-}(X^{i_{1},i_{2}}_{t_{l}},t_{l})\phi^{\prime\prime}(X^{i_{1},i_{2}}_{t_{l}}/\varepsilon)\right]$
$\displaystyle+\frac{\nu}{\varepsilon^{2}}\sum_{(i_{1},i_{2})\in
D_{b}^{+}}h_{1}h_{2}\sum_{l=0}^{k}h\mathbb{E}\left[1_{\\{t_{k}-t_{l}<\zeta(X^{i_{1},i_{2}}\circ\tau_{t_{k}})\\}}K_{\delta}(x,X^{i_{1},i_{2}}_{t_{k}})\theta_{+}(X^{i_{1},i_{2}}_{t_{l}},t_{l})\phi^{\prime\prime}(X^{i_{1},i_{2}}_{t_{l}}/\varepsilon)\right],$
(5.8)

for $x\in D$, and $\hat{u}=0$ otherwise. Notice that in the formula above,
sums over $(i_{1},i_{2})\in D$ denote both sums over the boundary layer
lattices $D_{b}^{\pm}$ and the outer layer lattice $D_{o}$. However, the last
two summations are over boundary lattices $D_{b}^{\pm}$ due to the remark we
made above regarding the support of the cutoff function. Thus,

$\displaystyle A_{i_{1},i_{2}}$ $\displaystyle=h_{1}h_{2}\text{ or
}h_{0}^{2},$ $\displaystyle\omega_{i_{1},i_{2}}$
$\displaystyle=\omega(x_{b\pm}^{i_{1}i_{2}},0)\text{ or
}\omega(x_{o}^{i_{1}i_{2}},0),$ $\displaystyle G_{i_{1},i_{2};t_{l}}$
$\displaystyle=G(x_{b\pm}^{i_{1}i_{2}},t_{l})\text{ or
}G(x_{o}^{i_{1}i_{2}},t_{l}),$ (5.9)

for boundary and outer layers respectively.

Since at every step we have to compute the boundary stress $\theta_{\pm}$, we
do this computation, according to (4.3) and (4.30), by applying the
corresponding derivatives to the formula (5.1). Indeed, this gives valid
approximations for the derivatives of the velocity $u$ as we work with the
mollified kernel $K_{\delta}(x,y)$. The expression for the derivatives is
given in terms of a similar to (5.1) formula where we replace the kernel
$K_{\delta}$ by the corresponding derivatives. We use the formulae for the
derivatives $\dfrac{\partial k^{-}}{\partial x}$ and $\dfrac{\partial
k^{+}}{\partial x}$ given in Lemma 3.2 which we similarly regularise using the
cutoff function $f_{\delta}$.

To handle the expectations in the above representation, we propose the
following numerical schemes.

Numerical scheme 1. In this numerical scheme, we omit the expectations in
(5.1) and run independent Brownian motions in (5.7). Therefore, we update the
diffusions $X_{t_{k}}^{i_{1},i_{2}}$, starting at $x^{i_{1}i_{2}}$ when $k=0$,
according to

$X^{i_{1},i_{2}}_{t_{k+1}}=X^{i_{1},i_{2}}_{t_{k}}+h\hat{u}(X^{i_{1},i_{2}}_{t_{k}},t_{k})+\sqrt{2\nu}(B_{t_{k+1}}^{i_{1},i_{2}}-B_{t_{k}}^{i_{1},i_{2}}),$
(5.10)

for $k\geq 0$, where $B^{i_{1},i_{2}}$ are independent Brownian motions. Here

$\displaystyle\hat{u}(x,t_{k+1})$ $\displaystyle=\sum_{(i_{1},i_{2})\in
D}A_{i_{1},i_{2}}\omega_{i_{1},i_{2}}K_{\delta}(x,X^{i_{1},i_{2}}_{t_{k}})1_{\\{t_{k}<\zeta(X^{i_{1},i_{2}}\circ\tau_{t_{k}})\\}}$
$\displaystyle+\sum_{(i_{1},i_{2})\in
D}A_{i_{1},i_{2}}hK_{\delta}(x,X^{i_{1},i_{2}}_{t_{k}})\sum_{l=0}^{k}G_{i_{1},i_{2};t_{l}}1_{\\{t_{k}-t_{l}<\zeta(X^{i_{1},i_{2}}\circ\tau_{t_{k}})\\}}$
$\displaystyle+\frac{\nu}{\varepsilon^{2}}\sum_{(i_{1},i_{2})\in
D_{b}^{-}}h_{1}h_{2}hK_{\delta}(x,X^{i_{1},i_{2}}_{t_{k}})\sum_{l=0}^{k}\theta_{-}(X^{i_{1},i_{2}}_{t_{l}},t_{l})\phi^{\prime\prime}(X^{i_{1},i_{2}}_{t_{l}}/\varepsilon)1_{\\{t_{k}-t_{l}<\zeta(X^{i_{1},i_{2}}\circ\tau_{t_{k}})\\}}$
$\displaystyle+\frac{\nu}{\varepsilon^{2}}\sum_{(i_{1},i_{2})\in
D_{b}^{+}}h_{1}h_{2}hK_{\delta}(x,X^{i_{1},i_{2}}_{t_{k}})\sum_{l=0}^{k}\theta_{+}(X^{i_{1},i_{2}}_{t_{l}},t_{l})\phi^{\prime\prime}(X^{i_{1},i_{2}}_{t_{l}}/\varepsilon)1_{\\{t_{k}-t_{l}<\zeta(X^{i_{1},i_{2}}\circ\tau_{t_{k}})\\}},$
(5.11)

for $x\in D$, and $\hat{u}=0$ otherwise with $A_{i_{1},i_{2}}$,
$\omega_{i_{1},i_{2}}$, $G_{i_{1},i_{2};t_{l}}$ given in (5.1).

Numerical scheme 2 In this scheme, we replace expectations in (5.1) by
averages with independent Brownian motions. Thus, we start the processes
$X_{t_{k}}^{m;i_{1},i_{2}}$ at $x^{i_{1}i_{2}}$ when $k=0$ and update them
according to

$X^{m;i_{1},i_{2}}_{t_{k+1}}=X^{m;i_{1},i_{2}}_{t_{k}}+h\hat{u}(X^{m;i_{1},i_{2}}_{t_{k}},t_{k})+\sqrt{2\nu}(B_{t_{k+1}}^{m}-B_{t_{k}}^{m}),$
(5.12)

for $k\geq 0$, where $B^{m}$ are independent Brownian motions for
$m=1,\ldots,N$, and

$\displaystyle\hat{u}(x,t_{k+1})$ $\displaystyle=\sum_{(i_{1},i_{2})\in
D}A_{i_{1},i_{2}}\omega_{i_{1},i_{2}}\frac{1}{N}\sum_{m=1}^{N}K_{\delta}(x,X^{m;i_{1},i_{2}}_{t_{k}})1_{\\{t_{k}<\zeta(X^{m;i_{1},i_{2}}\circ\tau_{t_{k}})\\}}$
$\displaystyle+\sum_{(i_{1},i_{2})\in
D}A_{i_{1},i_{2}}h\frac{1}{N}\sum_{m=1}^{N}K_{\delta}(x,X^{m;i_{1},i_{2}}_{t_{k}})\sum_{l=0}^{k}G_{i_{1},i_{2};t_{l}}1_{\\{t_{k}-t_{l}<\zeta(X^{m;i_{1},i_{2}}\circ\tau_{t_{k}})\\}}$
$\displaystyle+\frac{\nu}{\varepsilon^{2}}\sum_{(i_{1},i_{2})\in
D_{b}^{-}}h_{1}h_{2}h\frac{1}{N}\sum_{m=1}^{N}K_{\delta}(x,X^{m;i_{1},i_{2}}_{t_{k}})\sum_{l=0}^{k}\theta_{-}(X^{m;i_{1},i_{2}}_{t_{l}},t_{l})\phi^{\prime\prime}(X^{m;i_{1},i_{2}}_{t_{l}}/\varepsilon)1_{\\{t_{k}-t_{l}<\zeta(X^{m;i_{1},i_{2}}\circ\tau_{t_{k}})\\}}$
$\displaystyle+\frac{\nu}{\varepsilon^{2}}\sum_{(i_{1},i_{2})\in
D_{b}^{+}}h_{1}h_{2}h\frac{1}{N}\sum_{m=1}^{N}K_{\delta}(x,X^{m;i_{1},i_{2}}_{t_{k}})\sum_{l=0}^{k}\theta_{+}(X^{m;i_{1},i_{2}}_{t_{l}},t_{l})\phi^{\prime\prime}(X^{m;i_{1},i_{2}}_{t_{l}}/\varepsilon)1_{\\{t_{k}-t_{l}<\zeta(X^{m;i_{1},i_{2}}\circ\tau_{t_{k}})\\}},$
(5.13)

for $x\in D$, and $\hat{u}=0$ otherwise.

Note that in described numerical schemes, the stopping times
$\zeta(X^{i_{1},i_{2}})$ denote the first boundary hitting times for processes
$X^{i_{1},i_{2}}$ (or $X^{m;i_{1},i_{2}}$ for the second scheme), while
$\zeta(X^{i_{1},i_{2}}\circ\tau_{t_{k}})$ are seen as the last hitting times.
We therefore check if $X^{i_{1},i_{2}}_{t_{k+1}}$ crosses the boundary
$\partial D$ at every step when updating it according to (5.10) and (5.12).
Thus, after the first boundary crossing by $X^{i_{1},i_{2}}$ we do not count
the corresponding term in the first sum in (5.1) and (5.1).

Notice also that in (5.1) and (5.1) we separate the sums that track the
history of processes $X^{i_{1},i_{2}}$, i.e.

$\sum_{l=0}^{k}G_{i_{1},i_{2};t_{l}}1_{\\{t_{k}-t_{l}<\zeta(X^{i_{1},i_{2}}\circ\tau_{t_{k}})\\}},$

and

$\sum_{l=0}^{k}\theta_{\pm}(X^{i_{1},i_{2}}_{t_{l}},t_{l})\phi^{\prime\prime}(X^{i_{1},i_{2}}_{t_{l}}/\varepsilon)1_{\\{t_{k}-t_{l}<\zeta(X^{i_{1},i_{2}}\circ\tau_{t_{k}})\\}},$

as well as the corresponding sums with $X^{m;i_{1},i_{2}}$ for the second
scheme. Thus, we do not have to compute these sums at every step, instead we
store the whole sums and update them at each step adding new summands. When we
have a boundary crossing by $X^{i_{1},i_{2}}$ or $X^{m;i_{1},i_{2}}$, we set
the stored sum to zero and start updating it anew.

###  5.2 Numerical experiments

In this subsection, we use the numerical schemes described above to conduct
numerical experiments for certain flows. We exploit the Numerical Scheme 1
which is summarised in Algorithm 1 below.

Algorithm 1 Random vortex algorithm.

Define: $\omega_{i_{1},i_{2}}$ and $G_{i_{1},i_{2};t_{l}}$

Initialise: $X^{i_{1},i_{2}}\leftarrow x^{i_{1}i_{2}}$

Initialise: $\text{ind}^{i_{1},i_{2}}\leftarrow 1$ $\triangleright$ indicators
of $X^{i_{1},i_{2}}$ hitting the boundary

Compute:
$S^{i_{1},i_{2}}\leftarrow\theta_{\pm}(X^{i_{1},i_{2}}_{b\pm},t_{0})\phi^{\prime\prime}(X^{i_{1},i_{2}}_{b\pm}/\varepsilon)$

for $k\leftarrow 0$ to $N$ do

Compute: $\Delta^{i_{1},i_{2}}\leftarrow
h\hat{u}(X^{i_{1},i_{2}},t_{k})+\sqrt{2\nu}(B_{t_{k+1}}^{i_{1},i_{2}}-B_{t_{k}}^{i_{1},i_{2}})$

Update: $X^{i_{1},i_{2}}\leftarrow
X^{i_{1},i_{2}}+\Delta^{i_{1},i_{2}}\cdot\text{ind}^{i_{1},i_{2}}$

if $X^{i_{1},i_{2}}$ crosses $\partial D$ then

Set: $X^{i_{1},i_{2}}\leftarrow$ the crossing point

Set: $\text{ind}^{i_{1},i_{2}}\leftarrow 0$

Set: $S^{i_{1},i_{2}}\leftarrow 0$

end if

Compute: $\theta_{\pm}(X^{i_{1},i_{2}}_{b\pm},t_{k})$

Update: $S^{i_{1},i_{2}}\leftarrow
S^{i_{1},i_{2}}+\theta_{\pm}(X^{i_{1},i_{2}}_{b\pm},t_{k})\phi^{\prime\prime}(X^{i_{1},i_{2}}_{b\pm}/\varepsilon)$

end for

Experiment 1. We consider the flow past the flat-plate, i.e.
$D=\mathbb{R}^{2}\setminus\\{(x_{1},0):x_{1}\geq 0\\}$. The lattice points are
given by (5.3) and (5.4) with $N_{0}=15,N_{1}=30,N_{2}=45$. Therefore, we use
961 points in the outer layer and 1,426 points for each of the two thin
boundary layers. With $H=6$ and $H_{0}=0.1$ determining the size of the domain
and the thickness of the thin layers respectively, we find the mesh constants
$h_{0}=\frac{H}{N_{0}}=0.4,h_{1}=\frac{H}{N_{1}}=0.2,h_{2}=\frac{H_{0}}{N_{2}}\approx
0.0022$. We set the parameter $\varepsilon=0.05$ for the boundary sums in this
simulation.

We initialise the vector field
$u_{0}(x_{1},x_{2})=(U_{0}|x_{2}|(6-|x_{2}|)/9,0)$ for $-H\leq x_{1},x_{2}\leq
H$, where $U_{0}=\frac{\nu Re}{L}=125$ for $\nu=0.1,L=6$ and $Re=7500$. Thus,
the initial vorticity
$\omega_{0}(x_{1},x_{2})=-U_{0}(-2|x_{2}|+6)/9\cdot\textrm{sgn}(x_{2})$, where
for $x_{2}=0$ and $x_{1}<0$ we assume $\omega_{0}(x_{1},x_{2})=0$, and the
initial boundary stress $\theta_{+}(x_{1},0)=-2U_{0}/3$ and
$\theta_{-}(x_{1},0)=2U_{0}/3$. The external force $G$ is taken to be
identically zero.

We conduct our simulation with time steps $h=0.01$, the results for several
times $t$ are shown in Figures 5.1, 5.2 and 5.3. In these figures, the
streamlines are coloured by the velocity magnitude while the background is
coloured by the vorticity value. We also plot the boundary stress $\theta_{+}$
and $\theta_{-}$ as functions of position at boundary in Figures 5.4 and 5.5.

Experiment 2. In this experiment, the initial velocity is again taken to be
$u_{0}(x_{1},x_{2})=(U_{0}|x_{2}|(6-|x_{2}|)/9,0)$ and $G\equiv 0$. The
lattice points and other parameters are taken as before, however, we alter the
representation for $\hat{u}$. We ignore the boundary sums in this case, that
is we compute

$\hat{u}(x,t_{k+1})=\sum_{(i_{1},i_{2})\in
D}A_{i_{1},i_{2}}\omega_{i_{1},i_{2}}K_{\delta}(x,X^{i_{1},i_{2}}_{t_{k}})1_{\\{t_{k}<\zeta(X^{i_{1},i_{2}}\circ\tau_{t_{k}})\\}}.$

Therefore, we have to use a variant of Algorithm 1 with the above
representation for $\hat{u}$. The results are shown in Figures 5.6, 5.7 and
5.8 for the outer flow and boundary flows respectively, and Figures 5.9 and
Figures 5.10 for boundary stress values.

We notice that this version of the numerical scheme still allows to capture
some boundary phenomena as the boundary flows display certain chaotic
behaviour. However, at large times $t$ one has large regions close to the
boundary where the velocity does not change direction. We conclude therefore
that for simulating the boundary flows it is essential to use the boundary
sums in the velocity representation as they contribute to the chaotic
turbulent motion.

Experiment 3. We conduct this experiment for the wedge domain
$D=\mathbb{R}^{2}\setminus\left\\{x:-\alpha\leq\arg x\leq\alpha\right\\}$ with
$\alpha=\frac{\pi}{4}$. The size of the domain and numbers of lattice points
are taken as before, i.e. $H=6,H_{0}=0.1$ and $N_{0}=15,N_{1}=30,N_{2}=45$.
Notice that in this case $h_{1}=\frac{\sqrt{2}H}{N_{1}}\approx 0.2828$ as the
length of the boundary components contained in the domain is $\sqrt{2}H$.

We take the initial velocity $u_{0}(x_{1},x_{2})=0$, and the external force
term $G(x_{1},x_{2})=2U_{0}\frac{x_{2}}{H^{2}}$ for $U_{0}=500$ independent of
time. For this simulation, we choose smaller time steps $h=0.001$. The results
of the simulation are shown in Figures 5.11, 5.12, 5.13 — notice that in the
figures, the upper and lower boundary flows are rotated by $-\frac{\pi}{4}$
and $\frac{\pi}{4}$ respectively. Also, see Figures 5.14, 5.15 for boundary
stress values.

## References

  * [1] Anderson, C. and Greengard, C. 1985. On vortex methods _. SIAM J. Numer. Anal_. $\mathbf{22}$ (3), 413-440.
  * [2] Balakumar, B.J. and Adrian R.J. 2007 Large- and very-large-scale motions in channel and boundary-layer flows. _Phil. Trans._ $\mathbf{365}$, 665–681.
  * [3] Chauhan, K. Philip, J., De Silva, C.C.M., Hutchins, N. and Marusic, I. 2014 The turbulent/non-turbulent interface and entrainment in a boundary layer. _J. Fluid Mech_. 742, 119–151.
  * [4] Chorin, A. J. 1973 Numerical study of slightly viscous flow. _J. Fluid Mech._ $\mathbf{57}$, 785-796.
  * [5] Constantin, P. 2001 An Eulerian-Lagrangian approach for incompressible fluids: local theory. _J. Amer. Math. Soc._ $\mathbf{14}$ no. 2, 263-278 (electronic).
  * [6] Constantin, P. 2001 An Eulerian-Lagrangian approach to the Navier-Stokes equations. _Comm. Math. Phys._ $\mathbf{216}$, no. 3, 663-686.
  * [7] Constantin, P. and Iyer, G. 2011 A stochastic-Lagrangian approach to the Navier-Stokes equations in domains with boundary, _Ann. Appl. Probab_. 21, 1466-1492 (2011).
  * [8] Drivas, T.D. and Eyink, G.L. 2017 A Lagrangian fluctuation-dissipation relation for scalar turbulence. Part I. Flows with no boundary walls. _Journal of Fluid Mechanics_ , Volume $\mathbf{829}$, 25 October 2017 , pp. 153 - 189 DOI: https://doi.org/10.1017/jfm.2017.567
  * [9] Drivas, T. D. and Eyink, G.L. 2017 A Lagrangian fluctuation-dissipation relation for scalar turbulence. Part II. Wall-bounded flows. _Journal of Fluid Mechanics_. 829, 236-279 (2017).
  * [10] Eyink, G., Gupta, A., and Zaki, T. 2020 Stochastic Lagrangian dynamics of vorticity. Part 1. General theory for viscous, incompressible fluids. _Journal of Fluid Mechanics_ , 901, A2. doi:10.1017/jfm.2020.491
  * [11] Eyink, G., Gupta, A., and Zaki, T. 2020 Stochastic Lagrangian dynamics of vorticity. Part 2. Application to near-wall channel-flow turbulence. _Journal of Fluid Mechanics_ , 901, A3. doi:10.1017/jfm.2020.492
  * [12] Falkovich, G., Gawędzki, K. and Vergassola, M. 2001 Particles and fields in fluid turbulence, _Rev. Mod. Phys._ $\mathbf{73}$ 913-975.
  * [13] Feynman, R. P. 1948 Space-time approach to non-relativistic quantum mechanics. _Rev. Mod. Phys._ Vol. $\mathbf{20}$, No. 2, 367-387.
  * [14] Dawson, S.T.M. and McKeon, B.J. 2019 On the shape of resolvent modes in wall-bounded turbulence. _J. Fluid Mech_. $\mathbf{877}$, 682–716.
  * [15] Deardorff, J. W. 1970 A numerical study of three-dimensional turbulent channel flow at large Reynolds numbers. _J. Fluid Mech._ , $\mathbf{41}$:453-80
  * [16] DeGraaff, D. B. & Eaton, J. K. 2000 Reynolds-number scaling of the flat-plate turbulent boundary layer. _J. Fluid Mech_. $\mathbf{422}$, 319–346.
  * [17] Erm, L. P. & Joubert, P. N. 1991 Low Reynolds number turbulent boundary layers. _J. Fluid Mech_. $\mathbf{230}$, 1–44.
  * [18] Erm, L. P., Smits, A. J. & Joubert, P. N. 1985 Low Reynolds number turbulent boundary layers on a smooth flat surface in a zero pressure gradient. In _Proceedings of fifth Symposium on Turbulent Shear Flows_ , Ithaca, NY.
  * [19] Fletcher, C. A. J. 1991 _Computational techniques for fluid dynamics_ , Vol. I and II, second edition. Springer-Verlag.
  * [20] Friedman, A. 1964 _Partial differential equations of parabolic type_. Prentice-Hall, Inc.
  * [21] Goodman, J. 1987 Convergence of the random vortex method. _Comm. Pure Appl. Math._ $\mathbf{40}$(2), 189-220.
  * [22] Heisel, M., Dasari, T., Liu, Y., Hong, J., Coletti, F. & Guala, M. 2018 The spatial structure of the logarithmic region in very-high-Reynolds-number rough wall turbulent boundary layers. _J. Fluid Mech_. $\boldsymbol{857}$, 704–747.
  * [23] Head, M. R. & Bandyopadhyay, P. 1981 New aspects of turbulent boundary layer structure. _J. Fluid Mech_. $\mathbf{107}$, 297–338.
  * [24] Honkan, A. and Andreopoulos, Y. 1997 Vorticity, strain-rate and dissipation characteristics in the near-wall region of turbulent boundary layers. _J. Fluid Mech_. $\mathbf{350}$, 29–96.
  * [25] Kac, M. 1949 On Distributions of Certain Wiener Functionals. _Transactions of the American Mathematical Society_ , Jan., 1949, Vol. $\mathbf{65}$, No. 1 (Jan., 1949), pp. 1-13
  * [26] Keller, H. B. 1978 Numerical methods in boundary-layer theory. _Ann. Rev. Fluid Mech_. $\mathbf{10}$, 417-33.
  * [27] Lesieur, M., Métais, O. and Comte, P. 2005 _Large-Eddy simulations of turbulence_. Cambridge University Press.
  * [28] Li, J., Qian, Z., Xu, M. 2023 Twin Brownian particle method for the study of Oberbeck-Boussinesq fluid flows. https://doi.org/10.48550/arXiv.2303.17260
  * [29] Lilly, D. K. 1967 The representation of small-scale turbulence in numerical simulation experiments. In H. H. Goldstine (Ed.), _Proc. IBM Scientific Computing Symp. on Environmental Sciences_ , pp. 195-210. Yorktown Heights, NY: IBM.
  * [30] Long, D. G. 1988 Convergence of the random vortex method in two dimensions. _J. of Amer. Math. Soc_. $\mathbf{1}$(4 ), 779-804.
  * [31] Majda, A. J. and Bertozzi A. L. 2002 _Vorticity and incompressible flow_. Cambridge University Press.
  * [32] Moin, P. and Mahesh, K. 1998 Direct numerical simulation: a tool in turbulence research. _Annu. Rev. Fluid Mech._ $\mathbf{30}$, 539–78
  * [33] Orszag S.A. and Patterson G.S. 1972 Numerical simulation of three-dimensional homogeneous isotropic turbulence. _Phys. Rev. Lett_. $\mathbf{28}$, 76–79.
  * [34] Pope, S. B. 2000 _Turbulent flows_. Cambridge University Press.
  * [35] Prandtl, L. 1904 Über Flüssigkeitsbewegung bei sehr kleiner Reibung. _Proc. Third Intern. Math. Congress_ , Heidelberg, 848-491.
  * [36] Qian, Z. 2022 Stochastic formulation of incompressible fluid flows in wall bounded regions. https://doi.org/10.48550/arXiv.2206.05198
  * [37] Qian, Z., Süli, E. and Zhang, Y. 2022 Random vortex dynamics via functional stochastic differential equations. _Proc. R. Soc. A_ 478: 20220030. https://doi.org/10.1098/rspa.2022.0030
  * [38] Qian, Z., Qiu, Y., Zhao, L. and Wu, J. 2022 Monte-Carlo simulations for wall-bounded fluid flows via random vortex method. https://doi.org/10.48550/arXiv.2208.13233
  * [39] Rai, M. M. & Moin, P. 1993 Direct numerical simulation of transition and turbulence in a spatially evolving boundary layer._J. Comput. Phys_. $\mathbf{109}$, 169–192.
  * [40] Schlichting, H. and Gersten, K. 2017 _Boundary-Layer Theory_ (Ninth Edition). Springer.
  * [41] Spalart, P. R. 1988 Direct simulation of a turbulent boundary layer up to $R_{\theta}=1410$. _J. Fluid Mech_. $\mathbf{187}$, 61–98.
  * [42] Spalart, P. R. and Watmuff, J. H. 1993 Experimental and numerical study of a turbulent boundary layer with pressure gradients. _J. Fluid Mech_. $\mathbf{249}$, 337–371.
  * [43] Taylor, G. I. 1921 Diffusion by continuous movements. _Proc. Lond. Math. Soc_. $\mathbf{20}$, 196.
  * [44] Weinan, E. and Liu, J.-G. 1996 Vorticity boundary condition and related issues for finite difference schemes. _J. of Comp. Phys_. $\mathbf{124}$, 368–382.
  * [45] Wesseling, P. 2001 _Principles of computational fluid dynamics._ Springer-Verlag Berlin Heidelberg.
  * [46] Wu, X., Jacobs, R., Hunt, J. C. R. and Durbin, P. A. 1999 Simulation of boundary layer transition induced by periodically passing wakes. _J. Fluid Mech_. $\mathbf{398}$, 109–153.
  * [47] Wu, X. and Moin, P. 2008 A direct numerical simulation study on the mean velocity characteristics in turbulent pipe flow. _J. Fluid Mech_. $\mathbf{608}$, 81–112.
  * [48] Wu, X. and Moin, P. 2009 Direct numerical simulation of turbulence in a nominally zero-pressure-gradient flat-plate boundary layer. _J. Fluid Mech_. $\mathbf{630}$, pp. 5–41.
  * [49] Wu, X., Moin, P. and Hickey, J. P. 2014 Boundary layer bypass transition. _Physics of Fluids_ , $\mathbf{26}$, 091104\.

(a) $t=0.01$

(b) $t=0.25$

(c) $t=0.5$

(d) $t=1.0$

Figure 5.1: The outer layer flow at different times $t$.

(a) $t=0.01$

(b) $t=0.25$

(c) $t=0.5$

(d) $t=1.0$

Figure 5.2: The upper boundary layer flow at different times $t$.

(a) $t=0.01$

(b) $t=0.25$

(c) $t=0.5$

(d) $t=1.0$

Figure 5.3: The lower boundary layer flow at different times $t$.

(a) $t=0.01$

(b) $t=0.25$

(c) $t=0.5$

(d) $t=1.0$

Figure 5.4: The stress applied to the lower boundary at different times $t$.

(a) $t=0.01$

(b) $t=0.25$

(c) $t=0.5$

(d) $t=1.0$

Figure 5.5: The stress applied to the upper boundary at different times $t$.

(a) $t=0.01$

(b) $t=0.25$

(c) $t=0.5$

(d) $t=1.0$

Figure 5.6: The outer layer flow at different times $t$.

(a) $t=0.01$

(b) $t=0.25$

(c) $t=0.5$

(d) $t=1.0$

Figure 5.7: The upper boundary layer flow at different times $t$.

(a) $t=0.01$

(b) $t=0.25$

(c) $t=0.5$

(d) $t=1.0$

Figure 5.8: The lower boundary layer flow at different times $t$.

(a) $t=0.01$

(b) $t=0.25$

(c) $t=0.5$

(d) $t=1.0$

Figure 5.9: The stress applied to the lower boundary at different times $t$.

(a) $t=0.01$

(b) $t=0.25$

(c) $t=0.5$

(d) $t=1.0$

Figure 5.10: The stress applied to the upper boundary at different times $t$.

(a) $t=0.001$

(b) $t=0.05$

(c) $t=0.1$

(d) $t=0.2$

Figure 5.11: The outer layer flow at different times $t$.

(a) $t=0.001$

(b) $t=0.05$

(c) $t=0.1$

(d) $t=0.2$

Figure 5.12: The upper boundary layer flow at different times $t$.

(a) $t=0.001$

(b) $t=0.05$

(c) $t=0.1$

(d) $t=0.2$

Figure 5.13: The lower boundary layer flow at different times $t$.

(a) $t=0.001$

(b) $t=0.05$

(c) $t=0.1$

(d) $t=0.2$

Figure 5.14: The stress applied to the lower boundary at different times $t$.

(a) $t=0.001$

(b) $t=0.05$

(c) $t=0.1$

(d) $t=0.2$

Figure 5.15: The stress applied to the upper boundary at different times $t$.