Random vortex dynamics via functional stochastic differential equations

The statement of the point vortex method may be traced back to the paper [50], while the construction of effective computational methods by using vortex motions for viscous fluid flows began with the important work [12] (cf. also [14] and [13]) in which several important ideas, which became the core of vortex methods, were introduced to turn point-vortex approximations into powerful simulation tools. While point-vortex approximations are not exactly equivalent to the underlying fluid dynamics equations, and their convergence to true solutions of the equations is hard to prove, a great deal of effort has been devoted to proving the convergence of various numerical schemes for the approximate solution of fluid flow equations based on vortex approximations.

The convergence of the random vortex method for two-dimensional (2D) incompressible viscous fluid flow was studied in [28], and the rate of convergence was improved in [36]. Both authors utilized the fact that for 2D fluid flows vortex motions have the special feature that vortices are transported without experiencing nonlinear stretching.

The study of 3D inviscid fluid flows based on vortex dynamics has been an important discipline; one can find a review of advances in this field in [39] and [18], for example. 3D vortex methods for inviscid fluid flows were first presented by [3], and the convergence of the method was subsequently also proved. In [6, 7], the convergence of 3D vortex methods based on a Lagrangian frame was proved, and a new vortex method without the calculation of the Lagrangian frame was later proposed in [4]. In [17], a vortex method without mollifying the singular Biot–Savart kernel was proposed. For viscous flows certain numerical methods for computing solutions to fluid dynamics equations based on vortex dynamics were studied, for example, in [33] and [48]. However, the exact random vortex dynamics equations for the 3D incompressible Navier–Stokes equations have remained elusive. There are various stochastic representations of solutions to the 3D Navier–Stokes equations in the existing literature. In [16], a stochastic Lagrangian representation for the velocity of an incompressible fluid flow is obtained by making use of stochastic flows and Leray’s projection. In [56] the existence and uniqueness of the solutions to the Navier–Stokes equation was studied by using stochastic flow methods, and a relationship between mild solutions and stochastic representations has been investigated in [43]. In [55], stochastic flows associated with the Jacobian of the Taylor diffusion are used to study the solutions of the Navier–Stokes equation, while this approach is intended to be used for numerical algorithms. The authors of [2, 10] proposed different approaches, and obtained McKean–Vlasov type stochastic differential equations, which however involve a time-reversal in the coefficients. This formulation is therefore not amenable for numerical simulations. Other probabilistic representations are possible; for example, in [45] solutions are expressed as functionals of a Markov process indexed by the nodes of a binary tree.

Over the last two decades several deterministic vortex methods were proposed to approximate the dynamics of viscous fluid flows. For example, a particle strength exchange method (PSE) was explored in [20] and [19], where the authors approximate the Laplace operator via an integral operator and solve the related finite-dimensional ordinary differential equations. The Diffusion Velocity Method (DVM) was proposed and studied by [24, 42] and [8], where an artificial velocity field was incorporated into the dynamics in order to handle the diffusion of vorticity. In [54] the dynamics of Brownian motion type particles moving in chaotic flows are studied. More recently, machine learning techniques were applied in models of fluid flows; see, for example, [34], [35] and [27]. These methods have been useful for solving the 3D Navier–Stokes equations; they have, however, not revealed the underlying dynamics of vortices.

In this paper we derive a system of random dynamical equations for the three-dimensional incompressible Navier–Stokes equations, which is an exact closure of Taylor’s diffusion equation, the equation of motion for the “imaginary” Brownian fluid particles. The novel closure equations consist of a stochastic differential equation involving the law of its solution, coupled with an ordinary functional differential equation. The closed Taylor diffusion equations may then be used to design a random vortex method for the numerical approximation of the velocity field in three-dimensional incompressible viscous fluid flows by using a Monte Carlo method.

It is of course natural to explore the possibility of extending the proposed approach to other important classes of fluid flows, such as wall-bounded flows and the dynamics of viscous fluid flows in boundary layers. At the time of writing, significant progress has been made in this direction. For example, we are able to establish, with some subtle modifications, a forward type Feynman–Kac formula, and thus obtain similar, though more complicated, vortex dynamics for wall-bounded flows. The most subtle aspect in relation to wall-bounded flows with the approach proposed herein is handling the shear stress at the wall. The corresponding results will be published in a separate paper. The method presented here may also be applied to study other nonlinear problems by using diffusions constrained to domains with boundaries (cf. [1] and [51], for example). The reader may wish to consult [25], [26] and the references therein for a similar study.

The paper is organised as follows. In Section 2 we describe the main ideas in our approach and establish the notation that will be used in subsequent sections. In Section 3 we introduce a useful tool, called the pinned diffusion measure, the conditional law of an irreversible diffusion given a terminal location, which enables us to perform integration-by-parts in the proof of our main result in Section 5. In Section 4 we develop another probabilistic tool, a forward type Feynman–Kac representation for solutions of the vorticity equation. In Section 5, by exploiting the tools developed in Sections 3 and 4, a closed and exact set of random vortex dynamics equations is obtained; its discretisation gives rise to a direct numerical scheme defined in Section 6. Numerical results and the assessment of the numerical method on model problems are discussed in Section 7.

The random vortex method for two-dimensional (2D) viscous incompressible fluid flows based on exact distributional stochastic differential equations has been studied in [36], whose approach makes use of a special feature of the vorticity equations in two space dimensions, that vortices in 2D fluid flows are simply transported without nonlinear stretching of vortex motions.

Let $u=(u^{1},u^{2},u^{3})$ be the velocity field of an incompressible fluid flow with viscosity $\nu>0$, and let $\omega:=\nabla\wedge u$ be its vorticity. Then, the evolution of $u$ is governed by the incompressible Navier–Stokes equations

and

The vorticity $\omega$ evolves according to the vorticity equation

where

In the equations (2.3), (2.4) it is the vorticity $\omega$ that is considered as the main fluid dynamic variable, and the velocity $u$ is determined by the second equation, (2.4), together with the incompressibility condition (2.2). The fundamental difference between 2D fluid flows and 3D fluids flows lies in the fact that for 2D flows (i.e., flows for which the third velocity component $u^{3}=0$, and the first two velocity components, $u^{1}$ and $u^{2}$, depend on the first two coordinates, $x^{1}$ and $x^{2}$, only) the nonlinear vorticity stretching term $(\omega\cdot\nabla)u$ appearing on the right-hand side of (2.3) vanishes identically; hence the vorticity $\omega=(0,0,\omega^{3})$, which is identified with the scalar function $\omega^{3}$, evolves according to the following parabolic equation:

This crucial fact allows one to close the stochastic differential equations of Taylor’s Brownian motion $X$ (cf. [53, 5, 6, 7] and [36]):

issued from an initial state $x$, where $B$ is a three-dimensional Brownian motion, to provide powerful random vortex methods for computing 2D fluid flows; cf. [18, 36] and [39] for details.

To describe our approach and the contributions in the present work, we shall recall the key steps in the derivation of the 2D random vortex method in more detail and establish the notation that we will use throughout the paper. Let us begin by recalling several important mathematical structures associated with the velocity vector field $u(x,t)$. We shall discuss these constructions for a general vector field $b(x,t)$ for reasons that will become clear in what follows.

For a time-dependent, bounded and Borel measurable vector field $b(x,t)$ (where $x\in\mathbb{R}^{d}$, $t\geq 0$, and for the physically most relevant case of interest to us here, the spatial dimension $d=3$), let

which is a one-parameter family of second-order elliptic differential operators on $\mathbb{R}^{d}$. While the assumption of boundedness of $b$ is technical and may be replaced by other technical (weaker) conditions, we will not pursue this issue further, as our key aim in this paper is to derive a distributional stochastic differential equation, which will be used to design random vortex methods. If there is no risk of confusion, $L_{b(x,t)}$ will be abbreviated as $L_{b}$, where the arguments of $b$ are suppressed. Unless otherwise stated, we shall assume that the drift vector field $b(x,t)$ is solenoidal (i.e., divergence-free), that is, for every $t$, the divergence $\nabla\cdot b(\cdot,t)$ vanishes identically (in the sense of distributions on $\mathbb{R}^{d}$). Under these assumptions, the formal adjoint operator $L_{b}^{\star}=L_{-b}$, and therefore $L_{b}^{\star}$ is an elliptic operator of the same type as $L_{b}$. Let $\varGamma_{b}(x,t;\xi,\tau)$ (for $0\leq\tau<t$ and $\xi$, $x\in\mathbb{R}^{d}$) denote the fundamental solution to the forward heat operator $L_{b}-\frac{\partial}{\partial t}$, and $\varGamma_{b}^{\star}(x,t;\xi,\tau)$ (for $0\leq t<\tau$, $x$, $\xi\in\mathbb{R}^{d}$) the fundamental solution to the backward heat operator $L_{b}^{\star}+\frac{\partial}{\partial t}$, respectively (cf. Chapter 1 in [23] for their definitions, constructions and basic properties). Since $b$ is bounded, both $\varGamma_{b}(x,t;\xi,\tau)$ and $\varGamma_{b}^{\star}(x,t;\xi,\tau)$ are positive and continuous in $x,\xi$ and $t>\tau\geq 0$ (cf. [23] and [41], for example), and $\varGamma_{b}(x,t;\xi,\tau)=\varGamma_{b}^{\star}(\xi,\tau;x,t)$ for $t>\tau\geq 0$ and $x,\xi\in\mathbb{R}^{d}$. There are several probabilistic structures associated with a vector field $b(x,t)$ too, and these constructions will play an important rôle in this paper. Let $p_{b}(\tau,\xi,t,x)$ denote the transition probability density function of the diffusion process with associated infinitesimal generator $L_{b}$ (cf. [30] and [52]) in the following sense. Suppose that $X(\xi,\tau,t)$ (where $\xi\in\mathbb{R}^{d}$ is the initial state at instance $\tau$, $\tau\leq t$) is a diffusion process with associated infinitesimal generator $L_{b}$, i.e., $X$ is a (weak) solution to the stochastic differential equation (SDE)

and

where $B$ is a $d$-dimensional Brownian motion on a probability space $(\varOmega,\mathcal{F},\mathbb{P})$; $X(\xi,\tau,t)$ may be called the Brownian fluid particle with velocity $b(x,t)$. In the case when $\nu=0$, $t\rightarrow X(\xi,\tau,t)$ is just the integral curve of the vector field $b(x,t)$ starting at $\xi\in\mathbb{R}^{d}$ at time $\tau$. Then,

for any $t>\tau\geq 0$. When $\tau=0$, we will simply write $X(\xi,t)$ for $X(\xi,0,t)$ unless otherwise indicated. Since $b(x,t)$ is divergence-free, $L_{b}^{\star}=L_{-b}$; therefore,

for $t>\tau\geq 0$ and $\xi,x\in\mathbb{R}^{d}$.

Up to now we have collected several constructions that are valid in any number of space dimensions; we are now in a position to formulate the random vortex dynamics equations in two dimensions.

For a 2D viscous flow with velocity $u(x,t)$, the vorticity $\omega=\nabla\wedge u$ solves the forward “linear” parabolic equation (2.5), which may be written as

subject to appropriate boundary conditions at infinity. Hence,

for $(x,t)\in\mathbb{R}^{2}\times[0,\infty)$, provided that $\omega_{0}(x)$ tends to zero sufficiently fast as $|x|\rightarrow\infty$. Thanks to (2.10) the equality (2.12) may be written as

where $p_{u}(0,\xi,t,x)$ is the probability density function of the Taylor diffusion $X(\xi,t)$ with velocity $u$, which started at $\xi\in\mathbb{R}^{2}$ at $t=0$. Now recall that the Taylor diffusion is the weak solution to the SDE

for $t\geq 0$. The idea in random vortex methods is to rewrite the SDE (2.14) as a closed stochastic differential equation depending only on the initial datum $\omega_{0}$ and the distribution of Taylor’s diffusion $X$. To achieve this goal, one exploits the relations $\omega=\nabla\wedge u$ and $\nabla\cdot u=0$, which yields the Poisson equation $\Delta u=-\nabla\wedge\omega$, which, in turn, leads to the Biot–Savart law. Under the assumption that both $u$ and $\omega$ decay to zero at infinity sufficiently fast, we have that

where $K(z):=\frac{z^{\bot}}{|z|^{2}}$ is the 2D Biot–Savart kernel. Together with (2.13) one may rewrite

Thanks to this representation for the velocity $u(x,t)$ one is therefore able to close the stochastic differential equation (2.14):

and

for $(\xi,t)\in\mathbb{R}^{2}\times[0,\infty)$. The system of (2.16), (2.17) is a stochastic differential equation dependent on the law of its solution, i.e., it is a system of distributional stochastic differential equations.

There are obvious advantages to using numerical methods based on the exact stochastic dynamical system (2.16), (2.17). By employing a discreterisation of the finite-dimensional integral appearing in (2.16) and by appealing to the strong law of large numbers for the expectation, the closed system of stochastic differential equations (2.16), (2.17) leads to random vortex numerical schemes based on Monte Carlo simulation techniques for 2D incompressible fluid flows.

It would be helpful to formulate an exact random vortex dynamics system along the same lines for viscous incompressible fluid flows in three space dimensions. There is, however, a substantial obstacle one has to overcome, caused by the appearance of the nonlinear stretching term in the 3D vorticity equation. The nonlinear stretching term plays a particularly important role in 3D; indeed, in its absence there is no turbulent motion. Let us explain the key difficulty we have to overcome.

We begin with the vorticity equations (2.3), (2.4) for incompressible fluid flow with viscosity $\nu>0$ and velocity $u(x,t)$, which are rewritten as follows:

where $S=(S_{j}^{i})$ is the symmetric velocity gradient tensor, with entries $S_{j}^{i}:=\frac{1}{2}\left(\frac{\partial u^{i}}{\partial x^{j}}+\frac{\partial u^{j}}{\partial x^{i}}\right)$, and $\omega:=\nabla\wedge u$ is now, unlike in the 2D setting, a three-component vector field. In three space dimensions, the Biot–Savart law may be written as

where $K=(K^{1},K^{2},K^{3})$ is the singular vector kernel with components

for $i=1,2,3$. Because of the nonlinear stretching term, there is now no explicit representation for the vorticity $\omega$, which is considered as a solution to (2.18), in terms of the initial vorticity $\omega_{0}$ and the fundamental solution of $L_{-u}$, although the following implicit representation holds:

this expression, however, is not useful for the purpose of closing the system of stochastic differential equations (2.14) for Taylor diffusion. New ideas are therefore needed in order to formulate an exact system of random vortex dynamics.

A natural approach is to consider the symmetric velocity gradient tensor $S=(S_{j}^{i})$ as a multiplicative factor and apply the Feynman–Kac formula to the initial-value problem (2.18) by constructing a gauge functional from $S$ (cf. [22, 31, 32] for the linear case and [46, 47] for semilinear parabolic equations). Since the elliptic operator $L_{-u(x,t)}$ is temporally nonconstant and irreversible (which is canonical for flows observed in nature), the Feynman–Kac formula for $\omega$ thus obtained involves the distribution of a backward Taylor diffusion. Therefore there are at least two substantial difficulties one needs to overcome. First, this backward type Feynman–Kac formula cannot be used to close the system of stochastic differential equations (2.14), as we have to run a time-reversed Taylor diffusion in order to formulate the Feynman–Kac formula. Second, there are currently no satisfactory numerical methods for the approximate solution of semilinear parabolic equations based on the nonlinear version of the (backward) Feynman–Kac formula. In fact, existing methods appeal to numerical solutions of the corresponding PDEs, which is then in conflict with our aim to compute numerical approximations of solutions to the Navier–Stokes equations via Monte Carlo simulations.

We overcome these difficulties by devising two powerful new tools, which provide the means to close the system of stochastic differential equations (2.14) for Taylor diffusion with velocity $u(x,t)$. Recall that in two space dimensions, in order to close the system of stochastic differential equations (2.14) by eliminating the velocity $u$, two key steps are involved. The first is to devise a representation of the vorticity $\omega$ by ‘solving’ the vorticity equation, and the second step is to apply the Biot–Savart law to obtain a representation of the velocity $u$ in terms of $\omega$, then perform integration by parts (in two space dimensions one has to do this at any fixed time $t$) relying on the fact that $u$ is divergence-free, whereby $L_{-u}=L_{u}^{\star}$). In our approach, in the three-dimensional case, we devise a new representation formula, called a forward Feynman–Kac formula, stated in (4.2) below, for the vorticity $\omega$ by using the Taylor diffusion rather than a time-reversed diffusion. We are able to achieve this by using the fact that $u(x,t)$ is divergence-free, which implies that there is a convenient duality between the $L_{u}$-diffusion and the $L_{-u}$-diffusion conditional on their initial and final states; cf. Theorem 1 below. We then close the system of stochastic differential equations (2.14) by performing integration by parts at the path space level, – a technique which also appears to be new.

In this section we introduce the first technical tool in our study, pinned diffusion measures and a time-reversal theorem. This is considered as a path-space version of the integration by parts technique. Let $\nu>0$ be a constant, let $b(x,t)$ be a Borel measurable and bounded vector field depending on a parameter $t\geq 0$, and let $L_{b}:=\nu\Delta+b\cdot\nabla$. Let $\varOmega:=C([0,\infty),\mathbb{R}^{d})$ be the path space of all continuous paths in $\mathbb{R}^{d}$. The canonical coordinate process on $\varOmega$ is the ordered family of coordinate mappings $w\in\varOmega$, with value $w(t)$ at instance $t\geq 0$, where $t$ runs through $[0,\infty)$. Let $\mathcal{F}_{t}^{0}:=\sigma\{w(s):s\leq t\}$ be the smallest $\sigma$-algebra relative to which all coordinate mappings $w\rightarrow w(s)$ are measurable for all $s\leq t$, and $\mathcal{F}^{0}:=\sigma\{w(s):s\geq 0\}$. We have thus turned the path space $C([0,\infty),\mathbb{R}^{d})$ into a measurable space $(\varOmega,\mathcal{F}^{0})$. By the diffusion process $X=(X_{t})_{t\geq 0}$, with associated infinitesimal generator $L_{b}$, we mean the Markov family of probability measures

where $\mathbb{P}_{b}^{\xi,\tau}$ is the unique probability measure on $(\varOmega,\mathcal{F}^{0})$ such that the finite-dimensional joint distribution

is given by

for all $k=0,1,2,\ldots$, and $\tau<t_{1}<\cdots<t_{k}$. Here, thanks to the relation (2.10), one may identify the transition probability density function $p_{b}(\tau,\xi,t,x)$ with the fundamental solution $\varGamma_{-b}(x,t;\xi,\tau)$, which is positive and jointly Hölder continuous.

If $\tau\geq 0$, $T>\tau$ are fixed, and $\xi,\eta\in\mathbb{R}^{d}$ are two given states, then $\mathbb{P}_{b}^{\xi,\tau\rightarrow\eta,T}$ denotes the distribution of the $L_{b}$-diffusion process $X=(X_{t})_{t\geq 0}$ conditional on the event that $X_{s}=\xi$ for $s\leq\tau$ and $X_{s}=\eta$ for $s\geq T$. $\mathbb{P}_{b}^{\xi,\tau\rightarrow\eta,T}$ is again a probability measure on the continuous path space $\varOmega$, called the pinned diffusion measure starting from $\xi$ at instance $\tau$ and ending at state $\eta$ at the terminal time $T$, which can be described as follows. It is known that the laws of a Markov process conditional on the terminal values at some fixed time $T>0$ are also Markovian; see, for example, Section 14 in [21] and [15]. According to eq. (14.1) in [21, pp. 161], we define

for $\tau<s<t<T$ and $x,y\in\mathbb{R}^{d}$. It is then easy to verify that

and, for $\tau<s<t<r<T$, one has that

Hence, $\left\{q_{b}(s,x,t,y):\tau\leq s<t<T\right\}$ defines a temporal inhomogeneous Markov family. The $L_{b}$-pinned diffusion measure or the $L_{b}$-diffusion bridge measure $\mathbb{P}_{b}^{\xi,\tau\rightarrow\eta,T}$ is by definition the conditional law of $\mathbb{P}_{b}^{\xi,\tau}\left[\,\cdot\,|w(T)=\eta\right]$. It turns out that $\mathbb{P}_{b}^{\xi,\tau\rightarrow\eta,T}$ is the unique probability measure on $(\varOmega,\mathcal{F}^{0})$ such that the finite-dimensional marginal law

is given by

where

is any finite partition of $[\tau,T]$, and

By definition, if $\tau<s<T$ and if $f$ is a bounded Borel measurable function $f$, we have that

and we may therefore conclude that

for $s\in[\tau,T)$.

If $\tau=0$ and $T>0$, then, via the canonical embedding, $\mathbb{P}_{b}^{\xi,\tau\rightarrow\eta,T}$ can be considered (by an abuse of notation) as the probability measure on $\varOmega$ of all continuous paths in $\mathbb{R}^{d}$ with running time from $0$ to $T$.

The duality between the laws of the $L_{b}$-diffusions and the laws of the $L_{-b}$-diffusions registered in the next theorem plays an important rôle in the proposed new 3D random vortex method.

The assertion of this theorem may be reformulated as follows: suppose that $(X_{t})$ is an $L_{b}$-diffusion and $(Y_{t})$ is an $L_{-b^{T}}$-diffusion on the probability space $(\varOmega,\mathcal{F},\mathbb{P})$; then,

for any bounded or positive Borel measurable function $F$ on $\varOmega$.

To prove this theorem, we need the following simple fact concerning fundamental solutions.

Proof of Theorem 1. By definition,

is equal to the measure

which may be written as the following:

By Lemma 2, the above finite-dimensional distribution can be written as

which is the same as

Now the conclusion follows, as a distribution is determined uniquely by its family of finite-dimensional marginal distributions. This completes the proof. $\Box$

In this section we develop our second key tool, a forward Feynman–Kac formula, which is the key in deriving the 3D random vortex dynamics equations. We shall continue to assume that $b(x,t)$ is a bounded, Borel measurable and divergence-free vector field, with $x\in\mathbb{R}^{d}$ and $t\geq 0$.

In this section we derive the 3D random vortex dynamics SDE. To this end, let $u(x,t)$ be the velocity field of a viscous incompressible fluid flow (with viscosity $\nu>0$) without an external force applied, which moves freely in $\mathbb{R}^{3}$, and let $\omega:=\nabla\wedge u$. We assume that the velocity $u(x,t)$ and the vorticity $\omega(x,t)$ decay sufficiently fast as $|x|\rightarrow\infty$ and are sufficiently regular. Recall the vorticity equation

where $S$ is the symmetric velocity gradient tensor associated with $u$. Since $u$ is divergence-free, it follows that $L_{-u}=L_{u}^{\star}$.