Analytical solutions and asymptotic behaviors to the vacuum free boundary problem for 2D Navier-Stokes equations with degenerate viscosity

1.   Introduction

The evolving boundary of a viscous fluid can be modeled by the following compressible Navier-Stokes free boundary problem:

Here, $\rho$, ${{\mathbf{u}}} = (u_{1}, u_{2})\in\mathbb{R}^{2}$ and $p = p(\rho)$ denote, respectively, the density, the velocity field, and pressure of the fluid, which are functions of the space and time variable $(\mathbf{x}, t)\in\mathbb{R} ^{2}\times\lbrack0, \infty)$; $\tilde{\Omega}(t)\subset\mathbb{R}^{2}$ represents the changing volume occupied by a fluid at time $t$. The model described in Eq (1.1) can be used to describe the boundary expansion of gaseous stars, liquid flow in pipes, atmospheric flow, ocean currents, air currents around aircraft, and so on. For the polytropic gases, the pressure satisfies the common $\gamma$-law hypothesis

where $K > 0$ is a fixed constant and $\gamma$ is the adiabatic exponent. The constant $\gamma = c_{p}/c_{v}$ is the ratio of the specific heats, where $c_{p}$, $c_{p}$ are the specific heats per unit mass under constant pressure and constant volume, respectively. The different values of $\gamma$ imply different physical significance ^[2], for example, $\gamma = 5/3$, $\gamma = 7/5$, and $\gamma\rightarrow1^{+}$ correspond to a monatomic gas, a diatomic gas, and heavier molecules, respectively. In particular, the fluid is called isothermal if $\gamma = 1$. In this paper, we assume viscosity tensor $\Psi$ in (1.1)₂ to be of the following form:

where $I_{2}$ is the $2\times2$ identity matrix, and for simplicity, we set the viscosity coefficients

where the constants $k_{i}$ ($i = 1, 2, 3)\ $satisfy that

Equation (1.1) is completed by the vacuum free boundary condition (or continuous density condition)

where $\Gamma(t)$ denotes the moving interface. In general, as in ^[13] by Guo and Xin, the viscosity tensor can usually be given by the following form:

where $\mu_{1}$ and $\mu_{2}$ are the Lamé viscosity coefficients, and the inequality (1.8) is derived from physical constraints. In fact, in order to study the expansion of the vacuum boundary by using the energy method and the Bresch-Desjardins equality

Guo and Xin ^[13] moreover assumed that

Therefore, the range of parameters of the viscosity coefficients in our assumptions (1.3)–(1.5) generalize the one in Eqs (1.7)–(1.10), and it is important that the Bresch-Desjardins equality condition (1.9) is not necessarily satisfied.

Due to its physical importance and computational complexity in physics and mathematics, the vacuum free boundary problems have been widely studied in recent years, and important progress about the local and global well-posedness theory of weak, strong, or classical solutions has been made for both inviscid and viscous flows; interested readers may refer to ^{[3,12,14,25,26,41]} and references therein. Encountering the strong degeneracy on the vacuum free boundary of the density, the usual method of hyperbolic equations cannot be applied directly. Therefore, it is still a challenging problem to obtain the global existence of the system. The local well-posedness was only established recently for compressible inviscid flows (cf. ^[4,5,18]) and for compressible viscous flows (cf. ^[6,17]). We mention that for the vacuum free boundary problem (1.1) in multidimensional space, most of the global existence results are related to spherically symmetric solutions (cf. ^{[11,15,16,22,24,39]}) or affine ones (cf. ^[30,33]).

On the other hand, as pointed out by Yuen ^[37], the construction of analytical or exact solutions is important in mathematical physics and applied mathematics, due to that it can further classify their nonlinear phenomena, and a lot of important progress has been made in recent years. To begin, for inviscid and non-rotational flows (i.e., $k_{1} = k_{2} = k_{3} = 0$), the radially symmetric solutions and related exact solutions for the Euler equations were established in ^[19,32] and references therein. For rotational flows in 2D space, Zhang and Zheng ^[40] constructed analytical solutions for the Euler equations with $\gamma = 2$, which were generalized by Yuen ^[37] to the case $\gamma > 1$. In 3D space, Yuen ^[36] also gives a class of exact, rotational, infinite energy solutions to compressible or incompressible Euler and Navier-Stokes equations, where the solutions are similar to the famous Arnold-Beltrami-Childress (ABC) flow. On the other hand, the blowup phenomena of solutions have also attracted many researchers' attention in recent years. Dong and Yuen ^[9] studied the blowup of radial solutions to the compressible Euler equations (with or without damping) on some fixed bounded domains by introducing some new averaged quantities. When considering the influence of self gravity, Makino ^[27] proved the blowup (core collapsing) solutions to the 3D Euler-Poisson equations for $\gamma = 4/3$. It was extended by Deng, Xiang, and Yang in ^[8] to the case $N\geq3$ and $\gamma = (2N-2)/N$, then was generalized by Yuen ^[35] to the case which allows viscosity or frictional damping. For more results on the Euler equations and related equations, one may refer to ^{[1,13,21,23,28,34,38]} and references therein.

For the Navier-Stokes equation (1.1) in $\mathbb{R}^{N}$ ($N\geq2$) with both vacuum free boundary and stress free conditions, Guo and Xin ^[13] constructed spherically symmetric analytical solutions when density-dependent viscosity coefficients satisfy $k_{1} = \gamma > 1$, $k_{2} = 0$, and $k_{3} = \gamma-1 > 0$; in particular, the large time expanding behaviors at an algebraic rate of the free boundary are tracked. It is seen that the role of rotation is unknown in ^[13], while the reference ^[37] explores the effect of rotation. Thus, in this paper, based on the results in ^[13,37], we choose $r = a\left(t\right)$ as the free boundary and construct a class of self-similar analytical solutions for the Navier-Stokes equation (1.1) in 2D space with more general viscosity coefficients satisfying (1.4) and (1.5), which allows the effect of rotation. Moreover, the large time behavior of the free boundary is shown to be linear with respect to time. These rotational solutions will provide some new information on our understanding of ocean vortices and reference examples for numerical methods.

2.   Materials and methods

The existence of a class of self-similar analytical solutions to the isentropic compressible Navier-Stokes equations with vacuum free boundary in polar coordinates is established. In particular, these results can apply to the 2D Navier-Stokes equations with constant viscosity and the Euler equations. The special exact solutions constructed in this paper could also be applied in simulations of computing flows.

Considering that in polar coordinates, the original free boundary problem is simplified, and the corresponding equations are spatially dependent only on the radial variable $r$ (see (3.7)), this allows the analytical solutions of the equations to be solved sequentially. First, a pair of solutions $(\rho, u^r)$ of Eq (3.7)₁ can be obtained by using the known conclusion of self-similar solutions. Second, the analytical expression of $u^{\phi}$ can be obtained by substituting $(\rho, u^r)$ into Eq (3.7)₃ of $u^{\phi}$. Finally, by subsuming $(\rho, u^r, u^ {\phi})$ into Eq (3.7)₂, we can obtain a second order nonlinear ODE (with parameters of rotational strength $\xi$, adiabatic exponent $\gamma$, and viscosity coefficients; see (3.12)) that the free boundary $a(t) = r$ should satisfy. Next, the key is to study the existence and asymptotic behavior of the solution of the second-order equation. The fixed point theorem and the standard continuation argument can be used to prove the global existence of the solution, and the asymptotic behavior depends on the intrinsic structure of the equation. Specifically, by constructing an appropriate function $h(t)$ (see (3.42)), the monotonically increasing property of $a(t)$ can be obtained after a certain time $t_0$; the large time asymptotic behavior of $a(t)$ and $a^{\prime}(t)$ can be proved by using contradiction and the convexity of the function.

In the following sections, we will first give the equivalent formulation of the original free boundary problem in polar coordinates and state the main theorems in Subsection 3.1, then prove the Theorems 3.1 and 3.2 in Subsections 3.2 and 3.3, respectively. Finally, in supplementary, we give an explicit expression for the viscous terms $\mathrm{div\Psi}$ in polar coordinates.

3.   Results

3.1.   Formulation in polar coordinates and main results

The circular fluid region $\Omega\left(t\right)\in\mathbb{R}^{2}$ surrounded by vacuum in polar coordinates can be described as

where $r = \sqrt{x_{1}^{2}+x_{2}^{2}}$, the center of the region (0, 0) is fixed, and the free boundary $r = a(t)$ satisfies

where the positive and bounded constant $a_{0}$ represents the initial location of the free boundary $a(t)$. The velocity field has the form in Eulerian coordinates:

or, equivalently, in polar coordinates:

where ${{\mathbf{e}}}_{r} = \frac{(x_{1}, x_{2})}{r}$ and ${{\mathbf{e}}}_{\phi} = \frac{(-x_{2}, x_{1})}{r}$ are the two orthogonal unit vectors along the radial and the angular directions, respectively. Hence, the dissipative term $\mathrm{div\Psi}$ in Eq (1.1)₂ in the polar coordinate system has the following form:

(In fact, a detailed derivation of Eq (3.5) is shown in the supplementary.) Thus, Eqs (1.1)–(1.6) can be rewritten in polar coordinates as follows:

or equivalently as

with the initial conditions

and the Dirichlet boundary condition on the center of the region and the vacuum boundary condition on the free boundary:

In the following, we will use $C$ to denote the universal positive constants, which only depend on $\gamma$, $k_{i}$ $\left(i = 1, 2, 3\right)$, and the initial data such as $a_{0}$, $a_{1}$, and $H_{0}$ appearing in Theorem 3.1, but are independent of $t$, and they may change from one line to another. The labels $''x\lesssim y''$ and $''x\sim y''$ represent $''x\leq Cy''$ and $C_{1}y\leq x\leq C_{2}y$, respectively. The main results read:

Theorem 3.1. The problem (3.7)–(3.9) has a global solution of the form

where constants $\gamma > 1$, $\widetilde{k} > 0$, and $\xi\in\mathbb{R}$ are two arbitrary constants, $k_1$–$k_3$ satisfies the condition (1.5), and the free boundary $a\left(t\right) \in C^{2}([0, +\infty))$ satisfies the following Emden equation:

with initial values

Remark 3.1. In Theorem 3.2 below (see (3.15)), we can see that $a(t)$ is strictly positive, so the expressions (3.10)–(3.12) are well-defined, although the function a(t) appears as the denominator therein. The two constants $a_{0}$ and $a_{1}$ in (3.13) represent the initial location and slope of $a(t)$. If one sets $r = 0$ in (3.10) with a fixed adiabatic index $\gamma$, then $\widetilde{k}$ can characterize the magnitude of the fluid center density, and $\xi$ in (3.11) can describe the magnitude of the rotation intensity.

Remark 3.2. In 3D space, Yuen ^[36] also gives a class of exact, rotational, infinite energy solutions to Euler equations for $\gamma > 1$ in the following form:

with $C$, $c_{0}$, $c_{1}$, and $c_{2}$ arbitrary constants. Comparing (3.14), (3.10), and (3.11), it is interesting to see that the density and velocity functions in (3.14) both grow to infinity as time approaches infinity if $c_{1} > 0$, while the ones in (3.10) and (3.11) both decay to zero. The difference may be caused by the fact that Yuen considers the analytical solution of the whole-space problem, while we consider a bounded region with vacuum free boundary.

Theorem 3.2. For the Emden equation (3.12) with the parameters constraint (1.5), it has a unique and positive solution $a(t)$ such that

where

$H_{0} = \frac {1}{2}\left(a_{1}^{2}+\xi^{2}a_{0}^{-2}+\frac{K\gamma\widetilde{k}}{\gamma-1} a_{0}^{-2\left(\gamma-1\right) }\right)$, and $\bar{C} = \max\left\{ a_{0}, (2H_{0})^{1/2}\right\}$. Furthermore, the large time behaviors of $a(t)$ and $a^{\prime}\left(t\right)$ can be described as follows:

with constant

Remark 3.3. The constant $C_{0}$ appears in (3.16) as well-defined by (3.56) and (3.57). Similar to the derivation of Eq (3.12), if two or three of the three viscosity coefficients are constants, the following two special Emden equations can be obtained:

Case (1): $\lambda_{1}\left(\rho\right) = k_{1}$, $\lambda _{2}\left(\rho\right) = k_{2}$, and $\lambda_{3}\left(\rho\right) = k_{3}\rho^{\gamma}$, then $a(t)$ satisfies that

Case (2): $\lambda_{i}\left(\rho\right) = k_{i}$ $(i = 1, 2, 3)$, then $a(t)$ satisfies that

By comparing Eqs (3.12), (3.18), and (3.19), it can be seen that viscosity does affect the structure of the Emden equation. Moreover, Theorems 3.1 and 3.2 also apply to Eqs (3.18) and (3.19), except that $u^{\phi}\left(r, t\right)$ in Eq (3.11) will be replaced by $u^{\phi}\left(r, t\right) = \frac{\xi}{a^{2}\left(t\right) }r$. We also remark that the initial-value problem of Navier-Stokes equations was studied in ^[10] ($k_{1} > 0$, $k_{2} = 0$, $k_{3} > 0$), where the Cartesian solutions of the system without symmetry in $\mathbb{R}^{N}$ $\left(N\geq1\right)$ives a special solution to the Eule are given there.

Remark 3.4. If one sets $k_{1} = \gamma > 1$, $k_{2} = 0$, and $k_{3} = \gamma-1 > 0$, then (3.12) reduces to

which can be seen as a generalization of Eq (40) studied in ^[13]for a spherically symmetric case with $\xi = 0$. Note that the Bresch-Desjardins equality (1.9) in energy estimate is important for the spherically symmetric case. Here, the global analytical solution can still be obtained in Theorem 3.1 by directly studying Eq (3.12) of a(t), even though the Bresch-Desjardins equality does not hold true.

Remark 3.5. Note that for the special solution in (3.10), the viscosity term $(u_{r}^{r}+\frac{u^{r}}{r})_{r} = (u_{r}^{\phi}+\frac{u^{\phi}} {r})_{r} = 0$ (see (3.7)_{2, 3}), thus (3.10) and $u^{r}\left(r, t\right) = \frac{a^{\prime }\left(t\right) }{a\left(t\right) }r$, $u^{\phi}\left(r, t\right) = \frac{\xi}{a^{2}\left(t\right) }r$ also gives a special solution to the Euler equations, where $a\left(t\right)$ satisfies (3.19). In fact, all these solutions belong to the affine solution, and the simplest affine solution (spherically symmetric) or the general affine one for isentropic/non-isentropic Euler equations have been established in ^[30,31,33]. The innovation here is that we obtain a class of affine solutions for viscous fluids (Navier-Stokes equations with variable viscosity coefficients) with the same property of linear growth of the vacuum boundary.

Remark 3.6. We also mention that the solution constructed in (3.10) and (3.11) satisfies the physical vacuum boundary conditions (see ^[20,29,38]). Indeed, it follows from (3.28), (3.29), (3.23), and (3.12) that

which gives that

Integrating the equation above with respect to the space variable $r$ over $(r, a(t))$ (with $0 < r < a(t)$) and using the vacuum boundary condition (3.9) yields that

This, together with (3.15), implies that the sound speed $c = \sqrt{p^{\prime}\left(\rho\right) }$ is $C^{1/2}$-Hölder continuous (with respect to $r$) across the vacuum boundary, which is called the physical vacuum boundary condition.

3.2.   Proof of Theorem 3.1

Now, we show the proof in polar coordinates by some direct calculations, i.e., we will show that, Eq (3.7) has a class of solutions in the following form:

with constants $\widetilde{k} > 0$, $\xi\in\mathbb{R}$, the radius $r\in\left[ 0, a\left(t\right) \right]$, $f\left(s\right)$ is an arbitrary $C^{1}$ function of self-similar variable $s = \frac{r}{a\left(t\right) }$, and positive $a(t)\in C^{2}$ satisfies the Emden equation (3.12).

To begin, one can substitute $\rho$ and $u^{r}$ in (3.22) and (3.23) into Eq (3.7)₁ to obtain

Next, inserting $\rho$, $u^{r}$, and $u^{\phi}$ in (3.22) and (3.23) into the left-hand side of Eq (3.7)₃, one has

In view of (3.22), the third term on the righthand side of the equation above can be rewritten as

In order to seek a solution $u^{\phi}$ satisfying that (3.24) = 0, similar to that in ^[7,13], we set

integrating it over $\left(s, 1\right)$ and using the boundary condition that $f\left(1\right) = 0$ (due to (3.9)) to get

Hence, (3.25) can be rewritten as follows:

Thus, inserting (3.27) and $b\left(t\right)$ in (3.23) into (3.24), one gets

which implies that (3.7)₃ holds. Finally, we substitute (3.22) and (3.23) into Eq (3.7)₂ to deduce that

which is equivalent to the following:

Obviously, Eq (3.29) is exactly the Emden equation (3.12) by taking (3.27) into account. So, (3.22) and (3.23) are solutions to system (3.7)–(3.9). The proof of Theorem 3.1 is complete.

3.3.   Proof of Theorem 3.2

Note that Eq (3.12) belongs to the following type of ODEs:

with $g\left(t\right) \in(0, 1]$, $g^{\prime}\left(t\right) \leq0$, and two constants $C_{1} > 0$, $C_{2} > 0$. Indeed, the corresponding items $g\left(t\right)$, $C_{1}$, and $C_{2}$ to Eq (3.12) are as follows:

3.3.1.   Existence of solutions to (3.30)

In this subsection, we will prove the global existence of solutions to Eq (3.30) by establishing the local existence and global a priori estimates using the standard continuity argument. To this end, one can rewrite (3.30) as follows:

then

which gives that

Notice the equivalence of (3.30) and (3.33). We have the following local existence lemma by using the contraction mapping principle as in ^[7,13]; thus, we omit the details here.

Lemma 3.1. (Local existence) For Eq (3.30) with $\gamma > 1$ and $\widetilde{k} > 0$, there exists a small $T$ such that (3.30) has a positive solution $a(t)$, which is unique in $C^{2}([0, T])$ and satisfies $0 < a_{0}/2\leq a(t)\leq2a_{0}$.

Lemma 3.2. (Global existence) The Emden equation (3.30) has a positive solution $a(t)$, which is unique in $C^{2}([0, +\infty))$ and satisfies (3.15):

where $\bar{C} = \max\left\{ a_{0}, (2H_{0})^{1/2}\right\}$, $\underline{a} = \max\left\{ \left(\frac{C_{1}}{2\left(\gamma-1\right) H_{0}}\right) ^{1/\left[ 2\left(\gamma-1\right) \right] }, \frac{|\xi|g^{1/2}\left(t\right) }{(2H_{0})^{1/2}}\right\}$, $C_{1}$ and $g\left(t\right)$ are given by (3.30), and $H_{0}$ is defined by (3.37).

Proof. Assume $a(t)\in C^{1}([0, T])$ is a solution to (3.30). We first prove the a priori estimate

Multiplying (3.30) by $a^{\prime}(t)$ yields

then it follows that

Now, we define $H\left(t\right)$ as follows:

which, together with (3.35), for all $t\in\lbrack0, T]$, gives that

where $H_{0} = \frac{1}{2}\left[ a_{1}^{2}+\xi^{2}a_{0}^{-2}+\frac{C_{1} }{\gamma-1}a_{0}^{-2\left(\gamma-1\right) }\right]$. Obviously, (3.36) and (3.37) imply that

Due to $a_{0} > 0$ and the continuity property, one derives from (3.38) that

Thus, (3.38) and (3.39) yield that

It follows that

where $\bar{C} = \max\left\{ a_{0}, (2H_{0})^{1/2}\right\}$. Thus, (3.34) follows from (3.40) and (3.41). Therefore, combining the local existence, the a priori estimates in (3.34), and the standard continuity argument, we know that Eq (3.30) has a globally defined positive solution $a(t)$ satisfying (3.15). Thus, the proof of Lemma 3.2 is complete.

3.3.2.   Monotonically increasing property of $a(t)$ after time $t_0$

Let us define

It follows from (3.32) and (3.15) that

and

According to the sign of initial value $h(0)$, there are roughly two kinds of profiles of $a(t)$.

If $h\left(0\right) < 0$, due to the monotonicity and continuity property of $h(t)$, (3.43) implies that $h(t)$ will increase in a time interval until some finite time $t_{0} > 0$ (If $t_{0} = +\infty$, (3.42) implies that

then it holds that

and, hence,

Insert (3.46) into (3.44) to get, for a suitably large $t^{\ast} > 0$, that

which contradicts with (3.45). So, $t_{0} < +\infty$ holds.) such that $h(t_{0}) = 0$, and $t_{0}$ can be determined by

Thus, after time $t_{0}$, (3.44) implies that $h(t)\geq h(t_{0})$, namely,

where $t_{0}$ is determined by (3.47).

If $h\left(0\right) \geq0$, it follows from (3.44) and (3.42) that

so $a(t)$ increases for all time. Thus, it follows from (3.48) and (3.49) that

3.3.3.   Asymptotic behaviors of $a(t)$ and $a^{\prime}(t)$

To begin, we derive from (3.50) and the monotone bounded principle that the limit $\underset{t\rightarrow +\infty}{\lim}a(t)$ exists and belongs to $\left[ \underline{a}, +\infty\right]$. Moreover, we can claim that

Otherwise, suppose that there holds

for a suitably large $t^{\ast} > 0$, then it follows from (3.33) and (3.52) that

which contradicts (3.40). So, the supposition (3.52) fails, and (3.51) is true.

Due to (3.51) and (3.40), the following fact holds:

and (3.30) gives that, for a suitably large $t_{1} > 0$,

which implies that $a(t)$ is convex in $(t_{1}, +\infty)$. Thus, (3.40), (3.50), (3.53), and the monotone bounded principle yield that

and it follows that, for a suitably large $t^{\ast} > 0$,

for some positive constant $C_{0}$ to be determined later. By (3.55) and (3.31), we know the following integrability:

Now, letting $t\rightarrow +\infty$ in (3.33) and noting (3.56), one gets that

as $t\rightarrow +\infty$. Thus, (3.16) and (3.17) follow from (3.55), (3.57), (3.31), and the the L'Hospital rule, and we finish the proof of Theorem 3.2.

4.   Conclusions and discussion

In this paper, we established a class of self-similar analytical solutions to the vacuum free boundary problem for 2D isentropic Navier-Stokes equations with degenerate viscosity and studied their linear growth asymptotic behaviors for a large time. Here are a few ideas that we think are worth investigating. First, the expression (3.11) in polar coordinates is equivalent to, in Eulerian coordinates,

Hence, $(u_1, u_2)$ in Eqs (4.1) and (4.2) belongs to the class of affine solutions (or vector solutions). Therefore, it is reasonable to guess that the results obtained in this paper can be generalized to the case of affine solutions without symmetry, which probably requires the use of matrix theory, curve integration, and other related theories, as has been done in ^[10,33]. Second, we have selected special viscosity coefficients that satisfy (1.4) and (1.5):

Thus, the question is whether it is possible to extend the range of parameters in Eq (4.3), or to investigate a more general form of the viscosity coefficient as follows:

with some constants $\theta_i > 0$ and $k_{i}$ $(i = 1, 2, 3)$. Finally, we point out that the ideas and methods used in this paper can also be used to study the analytic solution and its large time behavior of the three-dimensional free boundary problem. In particular, it could be of great interest to consider the three-dimensional formulation of the problem of a spherically symmetric expansion of a compressible medium in a vacuum, and these issues will motivate our future work.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was partially supported by Natural Science Research of Jiangsu Higher Education Institutions of China (Natural Science Foundation of Colleges and Universities in Jiangsu Province) (22KJB110011), and Doctoral Research Fund of Huaiyin Normal University (31LKQ00). The author would like to thank Professor Zhengguang Guo for helpful discussions.

Conflict of interest

The author declares there is no conflict of interest in relation to this article.

Supplementary

Noting the definition of viscous stress tensor $\Psi$ given by (1.3), its divergence can be calculated as follows:

Now, we calculate the terms in Eq (2) in the following three cases in the polar coordinate system.

(1) $\nabla\rho$ and $\operatorname{div}{\mathbf{u}}$

For a scalar functions $f\left(r, t\right) = f\left(x_{1}, x_{2}, t\right)$ with $r = \sqrt{x_{1}^{2}+x_{2}^{2}}$, the chain rule gives us that

then

Let the function $f$ in (4) be the density or pressure of the fluid. One will have

For the vector function velocity field ${\mathbf{u}}$ (see (3.3)), we can deduce that

(2) $\bigtriangleup{\mathbf{u}} = div(\nabla{\mathbf{u}})$ and $div(\nabla{\mathbf{u}}^{T})$

If we set

then it follows that

and, similarly,

So, by noting the definitions of ${\mathbf{e}}_{r}$ and ${\mathbf{e}}_{\phi}$ in (3.4), we have

Similarly, we derive from (7) that

and, thus,

(3) $\left(\nabla\rho\cdot\nabla\right) {\mathbf{u}}$

By direct calculations, one has

Substituting expressions (5)–(10) into (1) and (2) produces

and

Finally, inserting (11)–(13) into (1) gives (3.5) directly.