Self-similar source-type solutions to the three-dimensional Navier-Stokes equations

We consider the Burgers equation [9]

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}=\nu\frac{\partial^{2}u}{\partial x^{2}},$$

(2)

which satisfies static scale-invariance under

$$x\to\lambda x,t\to\lambda^{2}t,u\to\lambda^{-1}u.$$

This means that if $u(x,t)$ is a solution, so is $u_{\lambda}(x,t)\equiv\lambda u(\lambda x,\lambda^{2}t),$ for any $\lambda(>0)$. It is readily checked that

$$\|u_{\lambda}\|_{L^{p}}=\lambda^{\frac{p-1}{p}}\|u\|_{L^{p}},$$

which shows that the $L^{1}$-norm is scale-invariant.

Let us clarify the two kinds of critical scale-invariance. Type 1 scale-invariance is achieved when we use a dependent variable whose physical dimension is the same as $\nu$. Type 1 is deterministic in nature where the additional term arising in the governing equations under dynamic scaling is minimised in number. Type 2 instead is statistical in nature where the additional terms under dynamic scaling are maximised in number so that a divergence form is completed and the dynamically-scaled equations have the Fokker-Planck operator as the linearisation. In the former the dependent variable has the same physical dimension as kinematic viscosity, whereas in the latter the argument of the Hopf characteristic functional (the independent variable) has the same physical dimension as the reciprocal of kinematic viscosity [15]. This approach provides a viewpoint from which the problem appears in the simplest possible form.

Critical scale-invariance of type 1 is achieved with the velocity potential $\phi,$ which is defined by $u=\partial_{x}\phi$. If $\phi(x,t)$ is a solution, so is $\phi(\lambda x,\lambda^{2}t).$ Under dynamic scaling for the velocity potential $\phi(x,t)=\Phi(\xi,\tau),$ $\xi=\frac{x}{\sqrt{2at}},\;\tau=\frac{1}{2a}\log t$ we have

$$\frac{\partial\Phi}{\partial\tau}+\frac{1}{2}\left(\frac{\partial\Phi}{\partial\xi}\right)^{2}=a\xi\frac{\partial\Phi}{\partial\xi}+\nu\frac{\partial^{2}\Phi}{\partial\xi^{2}},$$

(3)

whose linearisation has the Ornstein-Uhlenbeck operator. This is called type 1 (deterministic) scale-invariance where the number of additional terms is minimised, that is, only the drift term remains. Under dynamic scaling for velocity $u(x,t)=\frac{1}{\sqrt{2at}}U(\xi,\tau)$, we find

$$\frac{\partial U}{\partial\tau}+U\frac{\partial U}{\partial\xi}=a\frac{\partial}{\partial\xi}\left(\xi U\right)+\nu\frac{\partial^{2}U}{\partial\xi^{2}},$$

(4)

whose linearisation is the Fokker-Planck equation. Here the zooming-in parameter $a(>0)$ has the same physical dimension as $\nu$ and is on the same order of it. With this type 2 (statistical) scale-invariance where the number of additional terms is maximised meaning that a divergence form is completed with the addition of $aU$ term. As it is a second-order equation it has two independent solutions, of which we will focus on the Gaussian one. See Appendix C for the other non-Gaussian kinds of solutions.

Equation (4) is exactly soluble and its steady solution is called the source-type solution [11, 12]:

$$U(\xi)=\frac{U(0)\displaystyle{\exp\left(-\frac{a\xi^{2}}{2\nu}\right)}}{1-\displaystyle{\frac{U(0)}{2\nu}\int_{0}^{\xi}}\exp\left(-\frac{a\eta^{2}}{2\nu}\right)d\eta}.$$

(5)

The name has come from the time zero asymptotics

$$\lim_{t\to 0}\frac{1}{\sqrt{2at}}U(\xi)=M\delta(x),$$

where $M\equiv\int_{\mathbb{R}^{1}}u_{0}(x)dx$ and $U(0)=\sqrt{\frac{8a\nu}{\pi}}\tanh\frac{M}{4\nu}\approx\sqrt{\frac{a}{2\pi\nu}}M\;\;(\mbox{for}\;M/\nu\ll 1).$ Observe that (5) is a near-identity transformation of the Gaussian function. See [10, 11, 12].

It is also known that for $u_{0}\in L^{1}$ we have

$$t^{\frac{1}{2}\left(1-\frac{1}{p}\right)}\left\|u(x,t)-\frac{1}{\sqrt{2at}}U(\xi)\right\|_{L^{p}}\to 0\;\;\mbox{as}\;\;t\to\infty,$$

where $1\leq p\leq\infty.$

The simplest method for solving (4) without linearisation is as follows. Rewrite the equation

	$\displaystyle\frac{U^{2}}{2}$	$\displaystyle=$	$\displaystyle a\xi U+\nu\frac{dU}{d\xi}$
		$\displaystyle=$	$\displaystyle\nu\exp\left(-\frac{a\xi^{2}}{2\nu}\right)\frac{d}{d\xi}\left(U\exp\left(\frac{a\xi^{2}}{2\nu}\right)\right).$

By changing variables to $\tilde{U}=U\exp\left(\frac{a\xi^{2}}{2\nu}\right),\eta=\frac{1}{2\nu}\int_{0}^{\xi}\exp\left(-\frac{a\zeta^{2}}{2\nu}\right)d\zeta,$ we find

$$\frac{d\tilde{U}}{d\eta}=\tilde{U}^{2},$$

which is readily integrable. Alternatively we may solve the equation (2.1) by regarding it as a Bernoulli equation.

It may be in order to comment on the significance of source-type solution. When we recast (5) as

$$U(\xi)=-2\nu\frac{\partial}{\partial\xi}\log\left({1-\displaystyle{\frac{U(0)}{2\nu}\int_{0}^{\xi}}\exp\left(-\frac{a\eta^{2}}{2\nu}\right)d\eta}\right),$$

(6)

which is reminiscent of the celebrated Cole-Hopf transform. In other words, the source-type solution encodes the vital information of the nonlinear term in the case of the Burgers equation. Note that the error-function itself $\int_{0}^{\xi}\exp\left(-\frac{a\eta^{2}}{2\nu}\right)d\eta$ in (5, 6) is a self-similar solution to the heat equation. This suggests that studying source-type solution of the Navier-Stokes equations may give a hint on how to characterise their long-time evolution by a heat flow.

The operator $L=\triangle^{*}\equiv\triangle+\frac{a}{\nu}\partial_{\xi}(\xi\cdot)$ is not self-adjoint. It is possible to find a function $G$¹¹1With a slight abuse of notation this $G$ should be distinguished from the Gaussian function used in Section 3. such that $L^{\dagger}G(\xi)=-\delta(\xi)$ holds, where $L^{\dagger}\equiv\partial_{\xi}^{2}-\frac{a}{\nu}\xi\partial_{\xi}$ is the adjoint of $L$. In fact $G(\xi)\propto D\left(\sqrt{\frac{a}{2\nu}}\xi\right),$ where $D(\cdot)$ denotes the Dawson’s integral, defined by $D(x)\equiv e^{-x^{2}}\int_{0}^{x}e^{y^{2}}dy$ [18]. However, because $G$ decays slowly at large distances $G(\xi)\propto 1/\xi$ as $|\xi|\to\infty$, it cannot be used as a Green’s function, at least, in the usual manner.

The inversion formula for $\triangle^{*}$ can be obtained by an alternative method. Recall that based on $\frac{1}{a}=\int_{0}^{\infty}e^{-at}dt\;(a>0),$ the fundamental solution to the Poisson equation in 1D is given by

$$(\nu\triangle)^{-1}\equiv-\int_{0}^{\infty}dse^{\nu s\triangle}=\frac{|\xi|}{2\nu}*,$$

where * denotes convolution. Likewise for the fundamental solution to the Fokker-Planck equation in 1D we write

$$(\nu\triangle^{*})^{-1}\equiv-\int_{0}^{\infty}dse^{\nu s\triangle^{*}}=\int_{-\infty}^{\infty}d\eta g(\xi,\eta),$$

where

$$g(\xi,\eta)\equiv\frac{-1}{\sqrt{2\pi\nu}}{\rm f.p.}\int_{\sqrt{a}}^{\infty}\frac{d\sigma}{\sigma^{2}-a}e^{-\frac{1}{2\nu}\left(\sigma\xi-\eta\sqrt{\sigma^{2}-a}\right)^{2}}$$

and f.p. denotes the finite part of Hadamard, e.g. [13, 14]. It can be verified by changing the variable from $s$ to $\sigma=\sqrt{\frac{a}{1-e^{-2a\tau}}},$ using the solution of the Fokker-Planck equation

$$e^{\nu\tau\triangle^{*}}f=\left(\frac{a}{2\pi\nu(1-e^{-2a\tau})}\right)^{1/2}\int_{\mathbb{R}^{1}}e^{a\tau}f(e^{a\tau}\eta)\exp\left(-\frac{a}{2\nu}\frac{(\xi-\eta)^{2}}{1-e^{-2a\tau}}\right)d\eta.$$

As we will consider the 3D Navier-Stokes equations, for which methods of exact solutions are unavailable, we treat (4) by approximate methods as an illustration. Because the inversion of $\triangle^{*}$ is unwieldy, we will seek a workaround by which we can dispense with it.

First we convert it to an integral equation by the Duhamel principle for the Fokker-Planck operator $\triangle^{*}$

	$\displaystyle U(\tau)=$	$\displaystyle e^{\nu\tau\triangle^{*}}U(0)-\int_{0}^{\tau}e^{\nu(\tau-s)\triangle^{*}}\partial\,\frac{U(s)^{2}}{2}ds$
	$\displaystyle=$	$\displaystyle e^{\nu\tau\triangle^{*}}U(0)-\int_{0}^{\tau}e^{\nu s\triangle^{*}}\partial\,\frac{U(\tau-s)^{2}}{2}ds.$

The long-time limit $U_{1}=\lim_{\tau\to\infty}e^{\nu s\triangle^{*}}U(0)$ is given by

$$U_{1}=\left(\frac{a}{2\pi\nu}\right)^{1/2}Me^{-\frac{a}{2\nu}\xi^{2}}\;\;\mbox{with}\;\;M=\int_{-\infty}^{\infty}U(0)d\xi.$$

We may consider a number of different iteration schemes. For example, the following option (1), also known as the Picard iteration, requires the inversion $(\triangle^{*})^{-1}$:

$$\mbox{Successive approximation (1):}\;\;U_{n+1}=U_{1}-\int_{0}^{\infty}e^{\nu s\triangle^{*}}\partial\,\frac{U_{n}^{2}}{2}ds,$$

$$\mbox{in particular, for}\;\;n=1:\;\;U_{2}=U_{1}-\int_{0}^{\infty}e^{\nu s\triangle^{*}}\partial\,\frac{U_{1}^{2}}{2}ds.$$

Note that $U_{n}$ is a steady function at each step.

Alternatively we first consider the steady equation

$$\triangle^{*}U\equiv\triangle U+\frac{a}{\nu}(\xi U)_{\xi}=\frac{1}{\nu}\left(\frac{U^{2}}{2}\right)_{\xi}$$

and then introduce iteration schemes:

$$\mbox{Iteration scheme (2a):}\;\;\triangle U_{n+1}+\frac{a}{\nu}(\xi U_{n+1})_{\xi}=\frac{1}{\nu}\left(\frac{U_{n}^{2}}{2}\right)_{\xi},\;\;(n\geq 0)$$

$$\mbox{For}\;\;n=1:\;\;\triangle U_{2}+\frac{a}{\nu}(\xi U_{2})_{\xi}=\frac{1}{\nu}\left(\frac{U_{1}^{2}}{2}\right)_{\xi},$$

or

$$\mbox{Iteration scheme (2b):}\;\;\triangle U_{n+1}=-\frac{a}{\nu}(\xi U_{n})_{\xi}+\frac{1}{\nu}\left(\frac{U_{n}^{2}}{2}\right)_{\xi},\;\;(n\geq 1)$$

$$\mbox{For}\;\;n=1:\;\;\triangle U_{2}=\underbrace{-\frac{a}{\nu}(\xi U_{1})_{\xi}}_{=\triangle U_{1}}+\frac{1}{\nu}\left(\frac{U_{1}^{2}}{2}\right)_{\xi}.$$

Note that iteration schemes (1) and (2a) coincide with each other at $n=1$.

For the Burgers equation we can work out the two kinds of approximations to the second-order analytically. After straightforward algebra they are

	$\displaystyle\mbox{(1)}\;\;U$	$\displaystyle\approx Ce^{-\frac{a}{2\nu}\xi^{2}}\left(1+\frac{C}{2\nu}\int_{0}^{\xi}e^{-\frac{a}{2\nu}\eta^{2}}d\eta\right),$
	$\displaystyle\mbox{(2b)}\;\;U$	$\displaystyle\approx Ce^{-\frac{a}{2\nu}\xi^{2}}+\frac{C^{2}}{2\nu}\int_{0}^{\xi}e^{-\frac{a}{\nu}\eta^{2}}d\eta,$

where $C\approx\sqrt{\frac{a}{2\pi\nu}}M.$ On this basis we estimate the strength of the nonlinear term $N(\xi).$ The source-type solution is a near identity transform of the Gaussian function. In its series expansion in the Reynolds number $Re=M/\nu$ after non-dimensionalisation, the nonlinear correction term has $Re$ as its prefactor. Consider the scheme (1), or (2a) equivalently, taking $a/2\nu=1$ without loss of generality. We have

$$U_{2}=\frac{M}{\sqrt{\pi}}e^{-\xi^{2}}\left(1+\frac{Re}{2\sqrt{\pi}}\int_{0}^{\xi}e^{-\eta^{2}}d\eta\right).$$

Separating out the $Re$-dependence, or equivalently assuming that $Re=1,$ we define $N(\xi)$ by

$$N(\xi)=\frac{1}{2\sqrt{\pi}}\int_{0}^{\xi}e^{-\eta^{2}}d\eta$$

so that $U_{2}\propto 1+ReN(\xi)$ holds. The strength of nonlinearity is given by

$$N=\sup_{\xi}N(\xi)=\frac{1}{4}.$$

The same goes for (2b) from the above expressions for the 1D Burgers equation. Altogether we find

$$\mbox{(1)}\;\;N=\frac{1}{4}=0.25,$$

$$\mbox{(2b)}\;N=\frac{1}{4\sqrt{2}}\approx 0.2.$$

We conclude that the typical strength of nonlinearity is $N=O(10^{-1})$ irrespective of the choice of schemes.

The source-type solution is basically a near-identity function of the Gaussian form. It has been seen how the source-type solutions show up in the long-time limit in one and two spatial dimensions in [15]. Here we will take a look at cases in three and higher dimensions. From the Cole-Hopf transform we have

$$U_{i}(\bm{\xi},\tau)=-2\nu\frac{\partial_{\xi_{i}}\int_{\mathbb{R}^{3}}\psi_{0}(\lambda\bm{\eta})\exp\left(-\dfrac{a}{2\nu}\dfrac{|\bm{\xi}-\bm{\eta}|^{2}}{1-e^{-2a\tau}}\right)d\bm{\eta}}{\int_{\mathbb{R}^{3}}\psi_{0}(\lambda\bm{\eta})\exp\left(-\dfrac{a}{2\nu}\dfrac{|\bm{\xi}-\bm{\eta}|^{2}}{1-e^{-2a\tau}}\right)d\bm{\eta}}.$$

As we are going for the type 2 scale-invariance, differentiating it twice we find

$$\partial_{\xi_{j}}\partial_{\xi_{k}}U_{i}(\bm{\xi},\tau)=-2\nu\frac{\partial_{\xi_{i}}\partial_{\xi_{j}}\partial_{\xi_{k}}\int_{\mathbb{R}^{3}}\psi_{0}(\lambda\bm{\eta})\exp\left(-\dfrac{a}{2\nu}\dfrac{|\bm{\xi}-\bm{\eta}|^{2}}{1-e^{-2a\tau}}\right)d\bm{\eta}}{\int_{\mathbb{R}^{3}}\psi_{0}(\lambda\bm{\eta})\exp\left(-\dfrac{a}{2\nu}\dfrac{|\bm{\xi}-\bm{\eta}|^{2}}{1-e^{-2a\tau}}\right)d\bm{\eta}}+\ldots$$

$$=-2\nu\frac{\lambda^{3}\int_{\mathbb{R}^{3}}\partial_{i}\partial_{j}\partial_{k}\psi_{0}(\lambda\bm{\eta})\exp\left(-\dfrac{a}{2\nu}\dfrac{|\bm{\xi}-\bm{\eta}|^{2}}{1-e^{-2a\tau}}\right)d\bm{\eta}}{\int_{\mathbb{R}^{3}}\psi_{0}(\lambda\bm{\eta})\exp\left(-\dfrac{a}{2\nu}\dfrac{|\bm{\xi}-\bm{\eta}|^{2}}{1-e^{-2a\tau}}\right)d\bm{\eta}}+\ldots$$

The denominator then tends to $K_{ijk}\exp\left(-\frac{a}{2\nu}|\bm{\xi}|^{2}\right)$ as $\tau\to\infty,$ where $K_{ijk}=\int_{\mathbb{R}^{3}}\partial_{i}\partial_{j}\partial_{k}\psi_{0}(\bm{\eta})d\bm{\eta},\;(i=1,2,3).$ Hence

$$\partial_{\xi_{j}}\partial_{\xi_{k}}U_{i}(\bm{\xi},\infty)=-2\nu\left(\frac{K_{ijk}\exp\left(-\frac{a}{2\nu}|\bm{\xi}|^{2}\right)}{F_{ijk}(\bm{\xi})}+\ldots\right),$$

where the function $F_{ijk}$ is to be determined such that $\partial_{i}\partial_{j}\partial_{k}F_{ijk}\propto\exp\left(-\frac{a}{2\nu}|\bm{\xi}|^{2}\right),\;\;{\color[rgb]{0,0,0}\mbox{(no summation implied)}}.$ We can thus take

$$F_{ijk}(\bm{\xi})=-\frac{K_{ijk}}{2\nu}\int_{0}^{\xi_{1}}\exp\left(-\frac{a\xi^{2}}{2\nu}\right)d\xi\int_{0}^{\xi_{2}}\exp\left(-\frac{a\eta^{2}}{2\nu}\right)d\eta\int_{0}^{\xi_{3}}\exp\left(-\frac{a\zeta^{2}}{2\nu}\right)d\zeta+1.$$

Therefore after collecting other terms of derivatives we find in three dimensions, say, with $(i,j,k)=(1,2,3),$

$$\frac{\partial^{2}U_{1}}{\partial\xi_{2}\partial\xi_{3}}=K_{123}\exp\left(-\frac{a}{2\nu}(\xi_{1}^{2}+\xi_{2}^{2}+\xi_{3}^{2})\right)\frac{1+R(\xi_{1},\xi_{2},\xi_{3})}{(1-R(\xi_{1},\xi_{2},\xi_{3}))^{3}},$$

(7)

where

$$R(\xi_{1},\xi_{2},\xi_{3})=\frac{K_{123}}{2\nu}\int_{0}^{\xi_{1}}\exp\left(-\frac{a\xi^{2}}{2\nu}\right)d\xi\int_{0}^{\xi_{2}}\exp\left(-\frac{a\eta^{2}}{2\nu}\right)d\eta\int_{0}^{\xi_{3}}\exp\left(-\frac{a\zeta^{2}}{2\nu}\right)d\zeta$$

denotes the Reynolds number. Because $R$ is small the expression (7) is near-Gaussian. It can be verified that $K_{123}=\sqrt{\dfrac{32a^{3}}{\pi^{3}\nu}}\tanh\dfrac{M_{123}}{16\nu},$ where $M_{123}=\int\frac{\partial^{2}U_{1}}{\partial\xi_{2}\partial\xi_{3}}d\bm{\xi}.$ We can also write

$$\frac{\partial^{2}U_{1}}{\partial\xi_{2}\partial\xi_{3}}=-2\nu\frac{\partial^{3}}{\partial\xi_{1}\partial\xi_{2}\partial\xi_{3}}\log(1-R(\xi_{1},\xi_{2},\xi_{3})),$$

which reflects the Cole-Hopf transform more directly, that is, $\bm{U}=\nabla_{\bm{\xi}}\phi,\;\phi=-2\nu\log(1-R(\bm{\xi})).$ See Appendix A for the general form in $n$-dimensions.

The Burgers vortex was originally introduced to represent the reaction of a vortex under the influence of the collective effect of surrounding vortices in the ambient medium. When we write the steady solution in velocity and vorticity using cylindrical coordinates

$$\bm{u}=(u_{r},u_{\theta},u_{\phi})=(-ar,v(r),2az),\;\bm{\omega}=(0,0,\omega(r)),$$

the solution takes the following forms

$$\omega(r)=\frac{a\Gamma}{2\pi\nu}\exp\left(-\frac{ar^{2}}{2\nu}\right),$$

$$v(r)=\frac{\Gamma}{2\pi r}\left(1-\exp\left(-\frac{ar^{2}}{2\nu}\right)\right),$$

where $\Gamma\equiv\int_{\mathbb{R}^{2}}\omega_{0}(\bm{x})d\bm{x}$ denotes the velocity circulation.

In two dimensions dynamic scaling transforms take the following form

$$\bm{\xi}=\frac{\bm{x}}{\sqrt{2at}},\;\;\tau=\frac{1}{2a}\log t,\bm{\psi}(\bm{x},t)=\bm{\Psi}(\bm{\xi},\tau),$$

$$\bm{u}(\bm{x},t)=\frac{1}{\sqrt{2at}}\bm{U}(\bm{\xi},\tau),\omega(\bm{x},t)=\frac{1}{2at}\Omega(\bm{\xi},\tau).$$

In the 2-dimensional case, in order to achieve the critical scale-invariance of type 2, we must choose the second spatial derivative of the stream function, which is the vorticity, as the unknown.

The scaled form of the vorticity equation in two dimensions reads

$$\dfrac{\partial\Omega}{\partial\tau}+\bm{U}\cdot\nabla\Omega=\nu\triangle\Omega+a\nabla\cdot(\bm{\xi}\Omega),$$

where $\Omega$ satisfies the type 2 scale-invariance. It is known that the self-similar solution under scaling has a mathematically identical form as the Burgers vortex above. Indeed in the scaled variables the above expression can be written

$$\Omega(\xi)=\frac{a\Gamma}{2\pi\nu}\exp\left(-\frac{a|\bm{\xi}|^{2}}{2\nu}\right),\,\bm{\xi}=\frac{\bm{x}}{\sqrt{2at}}.$$

Note that $\frac{1}{2at}\Omega(\xi)=\frac{\Gamma}{4\pi\nu t}\exp\left(-\frac{|\bm{x}|^{2}}{4\nu t}\right)$ is an exact self-similar decaying solution²²2When $a=0$ an exact decaying solution is obtained by formally replacing $a\to\frac{1}{2t},$ which is known as the Lamb-Oseen vortex. with the following property

$$\lim_{t\to 0}\Omega(\cdot)=\Gamma\delta(\bm{x}).$$

It also satisfies the following asymptotic property, for $\omega(\cdot,0)\in L^{1},$

$$t^{1-\frac{1}{p}}\left\|\omega(\bm{x},t)-\frac{1}{2at}\Omega(\bm{\xi})\right\|_{L^{p}}\to 0\;\;\mbox{as}\;\;t\to\infty,$$

where $1\leq p\leq\infty,$ see e.g. [16].

Because the source-type solution coincides with the linearised solution, the inhomogeneous terms on the right-hand side of the approximations at each order vanish identically. Hence there is no way to set up successive approximations that can capture non-zero nonlinear corrections. The strength of nonlinearity is identically zero; $N=0$.

We consider the 3D Navier-Stokes equations written in four different dependent variables. Starting from the vector potential and taking a curl successively $\bm{u}=\nabla\times\bm{\psi},\;\bm{\omega}=\nabla\times\bm{u},\;\bm{\chi}=\nabla\times\bm{\omega},$ we have

$$\left\{\begin{array}[]{l}\dfrac{\partial\bm{\psi}}{\partial t}=\dfrac{3}{4\pi}{\rm p.v.}\displaystyle{\int_{\mathbb{R}^{3}}}\dfrac{\bm{r}\times(\nabla\times\bm{\psi}(\bm{y}))\,\bm{r}\cdot(\nabla\times\bm{\psi}(\bm{y}))}{|\bm{r}|^{5}}\;{\rm d}\bm{y}+\nu\triangle\bm{\psi},\\ \vskip 5.69046pt\cr\dfrac{\partial\bm{u}}{\partial t}+\bm{u}\cdot\nabla\bm{u}=-\nabla p+\nu\triangle\bm{u},\\ \vskip 5.69046pt\cr\dfrac{\partial\bm{\omega}}{\partial t}+\bm{u}\cdot\nabla\bm{\omega}=\bm{\omega}\cdot\nabla\bm{u}+\nu\triangle\bm{\omega},\\ \vskip 5.69046pt\cr{\color[rgb]{0,0,0}\dfrac{\partial\bm{\chi}}{\partial t}=\triangle(\bm{u}\cdot\nabla\bm{u}+\nabla p)+\nu\triangle\bm{\chi}}\end{array}\right.$$

(8)

where $\bm{r}\equiv\bm{x}-\bm{y}$ and p.v. denotes a principal-value integral. We also have $\bm{\omega}=-\triangle\bm{\psi},\bm{\chi}=-\triangle\bm{u},$ because of the incompressibility condition. The derivation of (8)₁ can be found in [17]. The final fourth equation (8)₄ is obtained by taking the Laplacian of the velocity equations (8)₂. For the $\bm{\chi}$ equations we may alternatively take a curl on the vorticity equations to obtain a form of equation different from the final line in (8), which is useful for handling inviscid fluids (details to be found in Appendix B). Under dynamic scaling

$$\bm{\xi}=\frac{\bm{x}}{\sqrt{2at}},\;\;\tau=\frac{1}{2a}\log t,\bm{\psi}(\bm{x},t)=\bm{\Psi}(\bm{\xi},\tau),$$

$$\bm{u}(\bm{x},t)=\frac{1}{(2at)^{1/2}}\bm{U}(\bm{\xi},\tau),\bm{\omega}(\bm{x},t)=\frac{1}{2at}\bm{\Omega}(\bm{\xi},\tau),\,\mbox{and}\,\bm{\chi}(\bm{x},t)=\frac{1}{(2at)^{3/2}}\bm{X}(\bm{\xi},\tau),$$

the 3D Navier-Stokes equations in four different unknowns are transformed respectively to

$$\left\{\begin{array}[]{l}\dfrac{\partial\bm{\Psi}}{\partial\tau}=\dfrac{3}{4\pi}{\rm p.v.}\displaystyle{\int_{\mathbb{R}^{3}}}\dfrac{\bm{\rho}\times(\nabla\times\bm{\Psi}(\bm{\eta}))\,\bm{\rho}\cdot(\nabla\times\bm{\Psi}(\bm{\eta}))}{|\bm{\rho}|^{5}}\;{\rm d}\bm{\eta}+\nu\triangle\bm{\Psi}+a(\bm{\xi}\cdot\nabla)\Psi,\\ \vskip 5.69046pt\cr\dfrac{\partial\bm{U}}{\partial\tau}+\bm{U}\cdot\nabla\bm{U}=-\nabla P+\nu\triangle\bm{U}+a(\bm{\xi}\cdot\nabla)\bm{U}+a\bm{U},\\ \vskip 5.69046pt\cr\dfrac{\partial\bm{\Omega}}{\partial\tau}+\bm{U}\cdot\nabla\bm{\Omega}=\bm{\Omega}\cdot\nabla\bm{U}+\nu\triangle\bm{\Omega}+a(\bm{\xi}\cdot\nabla)\bm{\Omega}+2a\bm{\Omega},\\ \vskip 5.69046pt\cr{\color[rgb]{0,0,0}\dfrac{\partial\bm{X}}{\partial\tau}=\triangle\left(\bm{U}\cdot\nabla\bm{U}+\nabla P\right)+\nu\triangle\bm{X}+a\nabla\cdot(\bm{\xi}\otimes\bm{X}),}\end{array}\right.$$

(9)

where $\bm{\rho}\equiv\bm{\xi}-\bm{\eta}.$

It is to be noted that the coefficient of the linear term increases in number with the increasing order of derivatives and for the variable $\bm{\chi}$ a divergence is completed in the convective term. Observe that type 1 scale-invariance is achieved with $\bm{\Psi}$ and type 2 scale-invariance with $\bm{X}$.

Using the Duhamel principle we convert the scaled Navier-Stokes equations (9)₄ into integral equations

$$\bm{X}(\bm{\xi},\tau)=e^{\nu\tau\triangle^{*}}\bm{X}_{0}(\bm{\xi})+\int_{0}^{\tau}e^{\nu s\triangle^{*}}\triangle\left(\bm{U}\cdot\nabla\bm{U}+\nabla P\right)(\bm{\xi},\tau-s)ds.$$

Here $\triangle^{*}\equiv\triangle+\frac{a}{\nu}\nabla\cdot(\bm{\xi}\otimes\cdot)$ and the action of whose exponential operator is given by

$$\exp(\nu\tau\triangle^{*})f(\cdot)=\left(\frac{a}{2\pi\nu(1-e^{-2a\tau})}\right)^{3/2}\int_{\mathbb{R}^{3}}e^{3a\tau}f(e^{a\tau}\bm{y})\exp\left(-\frac{a}{2\nu}\frac{|\bm{\xi}-\bm{y}|^{2}}{1-e^{-2a\tau}}\right)d\bm{y}$$

(10)

for any function $f,$ as can be verified by combining the heat kernel and dynamic scaling transforms.

The inverse operator associated with the fundamental solution to the Fokker-Planck equation in 3D is defined by

$$(\nu\triangle^{*})^{-1}\equiv-\int_{0}^{\infty}dse^{\nu s\triangle^{*}}=\int d{\color[rgb]{0,0,0}\bm{\eta}}g(\bm{\xi},\bm{\eta}),$$

where the Green’s function is given by

$$g(\bm{\xi},\bm{\eta})\equiv\frac{-1}{(2\pi\nu)^{3/2}}{\rm f.p.}\int_{\sqrt{a}}^{\infty}\frac{\sigma^{2}d\sigma}{\sigma^{2}-a}e^{-\frac{1}{2\nu}|\sigma\bm{\xi}-\bm{\eta}\sqrt{\sigma^{2}-a}|^{2}}.$$

This can be verified by changing the variable from $s$ to $\sigma=\sqrt{\frac{a}{1-e^{-2a\tau}}}$ in the solution (10) to the Fokker-Planck equation.

We consider the steady solution $\bm{X}(\bm{\xi})$ in the long-time limit of $\tau\to\infty$

$$\bm{X}(\bm{\xi})=\bm{X}_{1}(\bm{\xi})+\lim_{\tau\to\infty}\int_{0}^{\tau}e^{\nu s\triangle^{*}}\triangle\left(\bm{U}\cdot\nabla\bm{U}+\nabla P\right)(\bm{\xi},\tau-s)ds,$$

where $\bm{X}_{1}=\mathbb{P}\bm{M}G$ denotes the leading-order approximation (to be made more explicit in next subsection). This is one form of integral equations we are supposed to handle.

On the other hand, steady equations are obtained by assuming $\partial/\partial\tau=0$ in (9)₄

$$\triangle^{*}\bm{X}\equiv\triangle\bm{X}+\frac{a}{\nu}\nabla\cdot(\bm{\xi}\otimes\bm{X})=-\frac{1}{\nu}\triangle\left(\bm{U}\cdot\nabla\bm{U}+\nabla P\right),$$

or

$$\bm{X}=-\frac{1}{\nu}\left(\bm{U}\cdot\nabla\bm{U}+\nabla P\right)-\frac{a}{\nu}\triangle^{-1}\nabla\cdot(\bm{\xi}\otimes\bm{X}).$$

This is yet another form of the steady Navier-Stokes equations after dynamic scaling. It is noted that one of the potential problems associated with the nonlinear term has been eliminated without having a recourse to the Green’s function. It is this virtually trivial fact that allows us to set up a simple successive approximation.

To summarise, the steady Navier-Stokes equations after dynamic scaling can be written as

$$\bm{X}=-\frac{1}{\nu}\mathbb{P}\left(\bm{U}\cdot\nabla\bm{U}\right)-\frac{a}{\nu}\triangle^{-1}\nabla\cdot(\bm{\xi}\otimes\bm{X}),$$

(11)

or, by $\bm{X}=-\triangle\bm{U},$ we can express it solely in terms of $\bm{X}$ as

$$\bm{X}=-\frac{1}{\nu}\mathbb{P}\left(\triangle^{-1}\bm{X}\cdot\nabla\triangle^{-1}\bm{X}\right)-\frac{a}{\nu}\triangle^{-1}\nabla\cdot(\bm{\xi}\otimes\bm{X}).$$

(12)

This is the set of equations that we need to solve.

In passing we note the following facts before proceeding to the specific results. By the definition of scaled variables it is easily seen that for $p\geq 1$

$$t^{\frac{3}{2}\left(1-\frac{1}{p}\right)}\left\|\bm{\chi}(\bm{x},t)-\frac{\bm{X}(\bm{\xi})}{(2at)^{3/2}}\right\|_{L^{p}}=\left\|\bm{X}(\bm{\xi},\tau)-\bm{X}(\bm{\xi})\right\|_{L^{p}}.$$

This means that if $\left\|\bm{X}(\bm{\xi},\tau)-\bm{X}(\bm{\xi})\right\|_{L^{p}}\to 0$ as $\tau\to\infty,$ we have

$$t^{\frac{3}{2}\left(1-\frac{1}{p}\right)}\left\|\bm{\chi}(\bm{x},t)-\frac{\bm{X}(\bm{\xi})}{(2at)^{3/2}}\right\|_{L^{p}}\to 0\;\;\mbox{as}\;\;t\to\infty.$$

That is about the long-time asymptotics. On the other hand as time-zero asymptotics we have

$$\frac{\bm{X}(\bm{\xi})}{(2at)^{3/2}}\to\mathbb{P}\bm{M}\delta\;\;\mbox{as}\;\;t\to 0,$$

where $\delta(\cdot)$ is the Dirac mass.