Study of a semi-open queueing network with hysteresis control of service regimes

1.   Introduction

Viscoelastic materials include a wide range of fluids with elastic properties, as well as solids with fluid properties. The models of viscoelastic fluids formulated by Oldroyd [51,52] (see [56,53,5] for alternative derivations and perspectives), in particular, the classical Oldroyd-B model, have been studied by many authors (see [10,79,69] and the references cited therein), since the pioneering work of Renardy [58] and Guillopé–Saut [15].

It is well-known that viscoelasticity is a material property that exhibits both viscous and elastic characteristics with deformation. In particular, an elastic fluid strains when it is stretched and quickly returns to its original state once the stress is removed. This means that the elasticity will have a stabilizing effect in the motion of viscoelastic fluids. Moreover, the larger elasticity is, the stronger this stabilizing effect will be. This stabilizing effect has been widely investigated by many authors. In particular, recently the stabilizing effect in viscoelastic fluids was mathematically verified by Jiang et.al. based on the following three-dimensional (3D) incompressible Oldroyd model, which includes a viscous stress component and a stress component for a neo-Hookean solid:

where the unknowns ${v}: = {v}(x,t)$, and ${U}: = {U}(x,t)$ denote the velocity, and deformation tensor (a 3$\times$3 matrix valued function) of the viscoelastic fluids, resp.. The three positive (physical) constants $\rho$, $\mu$ and $\kappa$ stand for the density, shear viscosity coefficient and elasticity coefficient, resp., where $\kappa$ can be defined by the ratio between the kinetic and elastic energies, see [48], and we call the term ${\kappa} \mathrm{div}(UU^{ \mathrm{T}})$ the elasticity. The superscript $\mathrm{T}$ denotes transposition. We call (1.1)$_1$ the momentum equations, and (1.1)$_2$ the deformation equations. The divergence-free condition (1.1)$_3$ represents incompressibility of the viscoelastic fluids.

For strong solutions of the both Cauchy and initial-boundary value problems for (1.1), the authors in [49,50,9,45] have established the global(-in-time) existence of solutions in various functional spaces whenever the initial data is a small perturbation around the rest state $(0, I)$, where $I$ represents an identity matrix. The global existence of weak solutions to (1.1) with small perturbations near the rest state is established by Hu–Lin [24]. It is still, however, a longstanding open problem whether a global solution of the equations of incompressible viscoelastic fluids exists for any general large initial data, even in the two-dimensional case. At present, many authors have also obtained a lot of mathematical results for the well-posedness problem of the corresponding compressible case, see [23,25,26,27] for examples.

In the investigation of the well-posedness problem of (1.1) with initial value condition and the following non-slip boundary value condition of velocity

where $\Omega\subset \mathbb{R}^3$ denotes a fluid domain, the following basic energy identity plays an important role:

However, it is not easy to see from the above basic energy identity how elasticity effects the motion of viscoelastic fluids. To clearly see the effect, we shall rewrite (1.1) in Lagrangian coordinates.

To this purpose, let the flow map $\zeta$ be the solution to

Here and in what follows, the notations $f^0$ or $f_0$ denotes the initial data of $f$, if $f$ is a function of time variable $t$. In Lagrangian coordinates, the deformation tensor $\tilde{U}(y,t)$ is defined by the Jacobi matrix of $\zeta(y,t)$:

When we study this deformation tensor in Eulerian coordinates, it is defined by

where we have assumed that $\zeta$ is invertible with respect to variable $y$ for any given $t$. Moreover, by virtue of the chain rule, it is easy to check that $U(x,t)$ automatically satisfies the deformation equation (1.1)$_2$. This means that the deformation tensor in Lagrangian coordinates can be directly represented by $\zeta$, if the initial data $U^0$ also satisfies

Next, we proceed to rewrite the elasticity in Lagrangian coordinates. For this purpose, we should define the matrix $\mathcal{A}: = (\mathcal{A}_{ij})_{3\times3}$ via $\mathcal{A}^{ \mathrm{T}} = (\nabla\zeta)^{-1}: = (\partial_{j}\zeta_i)_{3\times3}^{-1}$, and the differential operators $\nabla_{\mathcal{A}}$, $\mathrm{div}_{\mathcal{A}}$ and $\Delta_{\mathcal{A}}$ are defined as follows, for $1\leqslant i\leqslant 3$,

for vector functions $w: = (w_{1},w_{2},w_{3})^{ \mathrm{T}}$ and $f_{i}: = (f_{i1},f_{i2},f_{i3})^{ \mathrm{T}}$, where we have used the Einstein summation convention over repeated indices, and $\partial_{k}: = \partial_{y_{k}}$. Let $J = \det \nabla \zeta$. Obviously, $\det U |_{x = \zeta} = J$.

It is well-known that

By (1.6) and the relation $\mathcal{A}^{ \mathrm{T}} \nabla\zeta = I$, we have

Now we assume $\det\nabla \zeta^0 = 1$, then $\det \nabla \zeta = 1$ due to the divergence-free condition (1.1)$_3$. For the incompressible case, (1.7) thus reduces to

Let $\eta = \zeta-y$ and

By virtue of (1.1)$_1$, (1.1)$_3$ and (1.2)$_1$, the system of equations for $(\eta,u,q)$ in Lagrangian coordinates reads as follows

which couples with the boundary-value condition:

It is easy to derive from (1.8) and (1.9) the following basic energy identity:

Let $(\eta^0,u^0)$ denote the initial data of $(\eta,u)$. Then integrating the above identity over $(0,t)$ yields

which implies that

In particular, for given the initial (perturbation) mechanical energy $I^0$, we see from the above identity that

Noting that $\eta$ physically represents the displacement function of particles around the rest state, hence the inequality (1.10) obviously exhibits that the elasticity has stabilizing effect in the flow.

Based on the above analysis result, Jiang et.al. mathematically proved that, under the proper large elasticity coefficient, the elasticity can inhibit flow instabilities such as Rayleigh–Taylor (RT) instability and thermal (or convective) instability. In addition, they further established the existence of large solutions of some class for the system of equations (1.1). We will review their results one by one.

2.   Viscoelastic Rayleigh–Taylor problem

2.1.   Stratified VRT problem in Eulerian coordinates

Consider two completely plane-parallel layers of stratified (immiscible) pure fluids, the heavier one is on top of the lighter one and both are subject to the earth's gravity. It is well-known that such an equilibrium state is unstable to sustain small disturbances, which will grow with time and lead to a release of potential energy, as the heavier fluid moves down under the gravitational force, and the lighter one is displaced upwards. This phenomenon was first studied by Rayleigh [57] and then Taylor [68] and is called therefore the RT instability. In the last decades, this phenomenon has been extensively investigated in mathematical, physical and numerical communities, see [7,71,17] for examples. It has been also widely investigated how other physical factors, such as rotation [7], internal surface tension [18,75,28], magnetic fields [35,29,30,34,73,74,32] and so on, influence the dynamics of the RT instability.

The RT instability was also often investigated in various models of viscoelastic fluids from the physical point of view, see [4] and [63] for examples. Later Jiang et.al. used the model (1.1) to investigate the effect of elasticity on the RT instability in viscoelastic fluids, see [36] and [39] for the stratified and non-homogeneous cases, resp.. Now we only introduce the RT instability in the stratified viscoelastic fluids.

To begin with, let us recall the RT problem of stratified viscoelastic fluids [36]:

Next, we should explain the notations in the above (stratified) VRT problem (2.1).

For each given $t>0$, $d: = d(x_{ \mathrm{h}},t):\mathbb{T}\mapsto (-l,h )$ is a height function of a point at the interface in stratified (viscoelastic) fluids; $\Sigma(t)$ is the set of interface defined by

where $l$, $h >0$,

$\mathbb{T}_i = 2\pi L_i(\mathbb{R}/\mathbb{Z})$, and $2\pi L_i$ ($i = 1$, $2$) are the periodicity lengths. Here and in what follows $f_{ \mathrm{h}}: = (f_1,f_2)$ for 3D vector $f$.

The subscripts $\pm$ in $f_\pm$ mean that a function or parameter $f_+$/$f_-$ is relevant to the upper/lower fluid. $\Omega_+(t)$ and $\Omega_-(t)$, defined by

and

are the domains of upper and lower fluids, resp.. $\Sigma_+$, resp. $\Sigma_-$ denote the fixed upper, resp. lower boundaries, are defined as follows:

For given $t>0$, $v_\pm(x,t): \Omega_\pm(t)\mapsto \mathbb{R}^3$, $U_\pm(x,t): \Omega_\pm(t)\mapsto \mathbb{R}^9$ and $p_\pm(x,t): \Omega_\pm(t)\mapsto \mathbb{R}$ are the velocities, the deformation tensors and the pressures of upper/lower fluids. $\rho_\pm$, $\kappa_\pm$ and $\mu_\pm$ are the density constants, viscosity coefficients and the elasticity coefficients of fluids, where the density of the upper fluid is heavier than the lower one, i.e.,

$g$ is the gravity constant. $\vartheta$ the surface tension coefficient. $\mathcal{S}_\pm$ represent the total strain tensors and are defined as follows:

where ${p}_\pm^g: = p_\pm + g\rho_\pm x_3$ and $\mathbb{D} v_\pm = \nabla v_\pm+\nabla v_\pm^{ \mathrm{T}}$.

For functions $f_\pm$ defined on $\Omega_\pm(t)$, we define $\left[\kern-0.15em\left[ { f_\pm } \right]\kern-0.15em\right] : = f_+|_{\Sigma(t)} -f_-|_{\Sigma(t)}$, where $f_\pm|_{\Sigma(t)}$ are the traces of the quantities $f_\pm$ on $\Sigma(t)$. $\nu$ is the unit outer normal vector on boundary $\Sigma(t)$ of $\Omega_-(t)$, and $\mathcal{C}$ twice of the mean curvature of the internal surface $\Sigma(t)$, i.e.,

Finally, we briefly explain the physical meaning of each identity in (2.1). The equations (2.1)$_1$–(2.1)$_3$ describe the motion of the upper heavier and lower lighter fluids driven by the gravitational field along the negative $x_3$-direction, which occupies the two time-dependent disjoint open subsets $\Omega_+(t)$ and $\Omega_-(t)$ at time $t$, resp.. The two fluids interact with each other by the motion of the free interface (2.1)$_4$ and the interfacial jump conditions (2.1)$_5$. The first jump condition in (2.1)$_5$ means that the velocity is continuous across the interface, while the second one is equivalent to

which represents that the jump in the normal stress is proportional to the mean curvature of the free surface multiplied by the normal to the surface [43,76]. The non-slip boundary condition of the velocities on both upper and lower fixed flat boundaries are described by (2.1)$_6$; and (2.1)$_7$–(2.1)$_8$ represent the initial status of the two fluids.

The problem (2.1) enjoys an equilibrium state (or rest) solution: $(v,U,{p}^g, d) = (0,I,\bar{p}^{g},\bar{d})$, where $\bar{d}\in (-l,h )$ is constant. We should point out that $\bar{p}^g$ depends on $x_3$ and can be uniquely (up a constant) evaluated by the relations

Without loss of generality, we assume that $\bar{d} = 0$ in this article. If $\bar{d}$ is not zero, we can adjust the $x_3$-coordinate to make $\bar{d} = 0$. Thus, $d$ can be regarded as the displacement away from the plane $\Sigma$.

To simply the representation of the problem (2.1), we introduce the indicator function $\chi_{\Omega_\pm(t)}$ and denote

To focus on the elasticity effect upon the RT instability, from now on we only consider the case $\vartheta = 0$, i.e., ignoring the effect of surface tension. Now, denoting the deviation of $(d,v,U,p^g)$ from the equilibrium state $(0,0,I,\bar{p}^g)$ by

one has a VRT problem in the perturbation form (with zero surface tension):

where $\mathcal{S}(\sigma,v,V +I)$ is defined by (2.3) with $(\sigma,v,V +I)$ in place of $(p ,v ,U )$, and we have denoted $\Omega(t): = \Omega_+(t)\cup\Omega_-(t)$, $\Sigma_-^+: = \Sigma_+\cup \Sigma_-$ and omitted the subscript $\pm$ in $\left[\kern-0.15em\left[ { f_\pm } \right]\kern-0.15em\right]$ for simplicity. Therefore, a zero solution is an equilibrium-state solution of the above perturbation VRT problem.

If ignoring the effect of elasticity (i.e., $\kappa = 0$) in (2.4) and deleting the deformation equation (2.4)$_2$, we obtain the classical RT problem of stratified immiscible pure fluids.

It is well-known that the classical RT problem is unstable, see [37,55] and the references cited therein. However, we will see that the VRT problem (2.4) may be stable due to the presence of elasticity.

2.2.   VRT problem in Lagrangian coordinates

It is well-known that the movement of the free interface $\Sigma(t)$ and the subsequent change of the domains $\Omega_\pm(t)$ in Eulerian coordinates will result in severe mathematical difficulties. There are three methods to circumvent such difficulties, i.e, the transform method of Lagrangian coordinates in [18], the method of translation transform with respect to $x_3$-variable in [55], and "geometric" reformulation method in [19]. Next, we shall modify and adopt the transform method of Lagrangian coordinates, so that the interface and the domains stay fixed in time. In particular, the VRT problem (2.1) in Lagrangian coordinates possesses a better mathematical structure.

Assume that there are invertible mappings

such that

and

where $\Sigma: = \mathbb{T} \times \{0\}$.

It should be noted that we define that $\Omega: = \Omega_+\cup \Omega_-$ in this section. We further define $\zeta: = \zeta_+\chi_{\Omega_+}+\zeta_-\chi_{\Omega_-}$, and the flow maps $\zeta$ as the solutions to

Denote the Eulerian coordinates by $(x,t)$ with $x = \zeta(y,t)$, whereas the fixed $(y,t)\in \Omega\times (0,\infty)$ stand for the Lagrangian coordinates. In addition, from incompressibility one gets $\det \nabla \zeta = 1$ in $\Omega$ as well as the initial condition (2.5).

Since $v_\pm$ and $\zeta_\pm^0$ are all continuous across $\Sigma$, one finds that $\Sigma(t) = \zeta_\pm(\Sigma,t)$. In view of the non-slip boundary value condition $v|_{\Sigma_-^+} = 0$, we have $y = \zeta(y, t)$ on $\Sigma_-^+$. We also assume that $\zeta (\cdot ,t)$ are invertible and $\Omega_{\pm}(t) = \zeta_{\pm}(\Omega_{\pm},t)$. Then we can define $\mathcal{A}$ as in Section 1 with $\zeta$ being provided by (2.6). We further define that

and refer to (1.4) for the definition of the operator $\nabla_{\mathcal{A}}$.

Now we further assume that $U^0_\pm = \nabla \zeta^0_\pm((\zeta^0_\pm)^{-1},0)$ in $\Omega_\pm(0)$. Setting $\eta = \zeta-y$ and the unknowns in Lagrangian coordinates

then, $(\eta,u,q)$ satisfies a so-called transformed VRT problem (without internal surface tension) in Lagrangian coordinates:

with $\tilde{U} (y,t): = \nabla \zeta (y,t)$. The operator $\mathrm{div}_{\mathcal{A}}$ is defined in (1.4).

2.3.   Main results

Before stating the main results for the transformed VRT problem (2.7), we shall introduce some notations of Sobolev spaces: for $k \geqslant 1$,

where $\Omega\!\!\!\!-\ : = \mathbb{T}\times (-l,h )$. Now we introduce the stability result for the transformed VRT problem (2.7).

Theorem 2.1. [36,Theorem 2.1] Under the stability condition

where

there is a sufficiently small constant $\delta >0$, such that for any $(u^0,\eta^0)\in H^{2}_0(\Omega)\times H^3_0(\Omega)$ satisfying

(1) the incompressible condition $\mathrm{div}_{\mathcal{A}^0}{u}^0 = 0$,

(2) the volume-preserving condition $\det(\nabla\eta^0+I) = 1$,

(3) $\|u^0\|_{H^2(\Omega)}+\| \eta^0 \|_{H^3(\Omega)} \leqslant\delta$,

(4) the initial data $(u^0, \eta^0)$ satisfies the compatibility condition

there exists a unique global solution $(u,\eta)\in C^0([0,\infty),H^{2}_0(\Omega)\times H^3_0(\Omega))$ to the transformed VRT problem (2.7) with an associated perturbation pressure $q$. Moreover, $(\eta,u,q)$ enjoys the following exponential stability estimate:

Here $\mathcal{A}^0: = (\nabla \eta^0+I)^{- \mathrm{T}}$, $\vec{n}^0: = {\mathcal{A}^0e_3}/{|\mathcal{A}^0e_3|}$, and the positive constants $\delta$, $\omega$ and $c$ depend on ${\rho}$, $\mu$, $g$, $\kappa$ and the domain $\Omega$.

There are some remarks of Theorem 2.1.

Remark 1. Theorem 2.1 still holds when $\Omega = \{x\in\mathbb{R}^3\; |\; -l<x_3<h \}$ (i.e., $L_1 = L_2 = +\infty$) or one of $\kappa_+$ and $\kappa_-$ is non-zero and the other is sufficiently large. In addition, if $\kappa_+ = \kappa_-$, then the stability condition is equivalent to $\kappa>\kappa_{ \mathrm{C}}$ with

Remark 2 The discriminant $\mathrm{C_r}$ enjoys the following estimate

where the positive constant $\tilde{c}$ is independent of $L_1$ and $L_2$. From the above estimate, we can conclude that

(1) $\mathrm{C_r}\to 0$ as $\kappa_+$ or $\kappa_-\to \infty$, which shows that a sufficiently large elasticity coefficient can prevent the RT instability from occurring. Here we present a reasonable physical explanation: when the equilibrium state is disturbed, the larger the elasticity coefficient is, the stronger the restoring force of the deformation tensor (i.e., elasticity force) is. Thus the sufficiently stronger elasticity can resist gravity and prevent the sagging of interface.

(2) $\mathrm{C_r}\to 0$ as $h$ or $l\to 0$, which also shows the heights $h$ and $l$ do influence the value of $\mathrm{C_r}$. In addition, the lengths $L_1$, $L_2$ do not predominantly affect the value of $\mathrm{C_r}$, when $g$, $\rho$, $l$ and $h$ are given. Compared with the classical RT problem of stratified immiscible fluids with internal surface tension, we know that the internal surface tension can constrain the growth of the RT instability, when the coefficient of surface tension $\vartheta$ is greater than a critical number $\vartheta_{ \mathrm{C}}: = g \left[\kern-0.15em\left[ { \rho} \right]\kern-0.15em\right] \max\{L^2_1,L^2_2\}$ [18,75], from which we see that the lengths $L_1$, $L_2$ do affect $\vartheta_{C}$, while the heights $h$, $l$ do not.

Remark 3. For the current viscoelastic RT problem, $\mathrm{C_r}>0$ due to $\left[\kern-0.15em\left[ { \rho} \right]\kern-0.15em\right] > 0$. However, if one further considers $\left[\kern-0.15em\left[ { \rho} \right]\kern-0.15em\right] \leqslant 0$, then $\mathrm{C_r}\leqslant 0$. In this case, Theorem 2.1 still holds. Moreover, Theorem 2.1 also holds for one-layer viscoelastic fluid in an infinite strip with the upper free boundary, which can be regarded as a spacial case of $\rho_+ = 0$ and $\rho_-\neq0$. This exponential stability result on a one-layer viscoelastic fluid also shows that the elasticity has stabilizing effect, because, for a one-layer viscous flow in an infinite strip with upper free boundary, Guo–Tice [20,19] used the two-tiers energy method to only show the existence of a unique solution that decays algebraically in time.

Remark 4. We remark on how to use the stability condition (2.8). Let us recall the basic energy identity of (2.7):

where $\tilde{\mathcal{A}}: = \mathcal{A}-I$, $\mathcal{E}_1 : = \|\sqrt{\rho}u\|^2_{L^2(\Omega)} -{E}_1 (\eta)$ and

We call $\mathcal{E}_1$ the (linear) mechanical energy and ${E}_1 (\eta)$ the (linear) potential energy.

Obviously, to derive the a priori stability (2.9), we shall first pose a stability condition, which makes sure the mechanical energy to be positive definite. It is easy to see that the condition (2.8) is just the stability condition, for which we look. In particular, under the stability condition, we have the following stabilizing estimate: there exists a positive constant $c$ depending on $g$, $\rho$ and $\kappa$ such that

Thanks to (2.10), we can use an energy method to establish Theorem 2.1.

The failure of the stability condition (2.8) results in the following instability result, which exhibits that the elasticity can not inhibit RT instability for small elasticity coefficient.

Theorem 2.2. [37,Theorem 1.1] Under the instability condition

a zero solution to the transformed VRT problem is unstable in the Hadamard sense, that is, there are positive constants $\Lambda$, $m_0$, $\epsilon$ and $\delta_0$, and a vector function $(\tilde{\eta}^0,\tilde{u}^0,\tilde{q}^0$, $\eta^ \mathrm{r},u^ \mathrm{r},q^ \mathrm{r})\in H^3_\sigma(\Omega)\times H^2_\sigma(\Omega)\times H^1(\Omega)\times H^3_0(\Omega)\times H_0^2(\Omega)\times H^1(\Omega)$, such that, for any $\delta\in (0,\delta_0)$ with

there is a unique strong solution $(\eta,u)\in C^0([0,T],H^{3}_{0}(\Omega)\times H^2_0(\Omega))$ to the transformed VRT problem with initial data $(\eta^0,u^0)$ and with a unique (up to a constant) associated pressure $q$ (with initial data $q^0$). In adiition, the solution satisfies

for some escape time $T^\delta: = \frac{1}{\Lambda} \mathrm{ln}\frac{2\epsilon}{m_0\delta}\in (0,T)$, where $\mathcal{A}: = (\nabla \zeta)^{-1}$, $\zeta: = \eta+y$, $\chi_{ \mathrm{h}}: = (\chi_1,\chi_2)^{ \mathrm{T}}$, $\chi = \eta$ or $u$ and $T$ denotes some time of existence of the solution $(\eta ,u)$. Moreover, the initial data $(\eta^0,u^0,q^0)$ satisfies $\det(\nabla \eta^0+I) = 1$ and necessary compatibility conditions:

where $\mathcal{A}^0: = (\nabla \eta^0+I)^{- \mathrm{T}}$ and $\vec{n}^0: = {\mathcal{A}^0e_3}/{|\mathcal{A}^0e_3|}$.

We give some remarks for Theorem 2.2.

Remark 5. In this review, we only introduce the results of stability and instability for the transformed VRT problem (2.7) without surface tension. Interested readers can further refer to [78,37] for the case with surface tension. In addition, the corresponding version of compressible case can be found in [31].

Remark 6. It should be noted that, if $\delta$ is sufficiently small, the solutions obtained in Theorems 2.1 and 2.2 automatically satisfy the following properties: for any given $t\geqslant 0$,

where $\zeta: = \eta+y$ and $\zeta_+(\Omega_+) \cup \zeta_-(\Omega_-)\cup \zeta_-(\Sigma) = \Omega\!\!\!\!\!-$. Exploiting the properties above, we can recover the exponential stability, resp. instability of the VRT problem (2.4) in Eulerian coordinates from Theorems 2.1, resp. 2.2 by an inverse transform of Lagrangian coordinates, see [37,Theorem 1.2] for an example.

Remark 7. Finally we remark on how to use the instability condition (2.11). The eigenvalue problem of the linearized VRT problem reads as follows:

We expect to look for the largest positive eigenvalue $\lambda>0$ and the corresponding eigenvalue function $(\tilde{\eta},\tilde{u}, \tilde{q})$. Since the above eigenvalue problem enjoys a fine symmetrical structure, the problem of looking for the largest positive eigenvalue reduces to the following variational problem:

where $\mathcal{A}: = \left\{w\in {H}_{\sigma }^1(\Omega)\; \left|\; \|\sqrt{\rho}w\|^2_{L^2(\Omega)} = 1\right.\right\}$. By (2.11), it holds that

This means that the variational problem (2.14) makes sense. Thus we can use the modified variational method in [18] to obtain the largest positive eigenvalue $\lambda$ and the corresponding eigenvalue function $(\tilde{\eta},\tilde{u})$ with an associated function $\tilde{q}$. In particular, we further define that

which is the solution with large exponent growth in time to the linearized transformed VRT problem. Thanks to this linear instability result, we can use a bootstrap instability method in [16] to establish Theorem 2.2.

It is well-known that, in the development of the RT instability, gravity first drives the third component $u_3$ of the velocity unstable, then the instability of the third component of the velocity further results in instability of the horizontal velocity $(u_1,u_2)$ and the height function $\eta_3(y_{ \mathrm{h}},0,t)$. Obviously, the instability results in (2.13) accords with this instability phenomenon. Moreover, the instability velocity is accompanied by the instability of displacement function of particles. However, elasticity can inhibit the growth of displacement function of particles. In particular, elasticity can recover the particles to equilibrium state under sufficiently large elasticity coefficient. This physical interpretation easily motivates to predict that the other flow instabilities, which cause the instability of displacement of particles, may be inhibited by elasticity. Next we further take the thermal instability in viscoelastic fluids as an example.

3.   Viscoelastic Rayleigh–Bénard problem

Thermal instability often arises when a fluid is heated from below. The classic example of this is a horizontal layer of fluid with its lower side hotter than its upper. The basic state is then one of rests with light and hot fluid below heavy and cool fluid. When the temperature difference across the layer is great enough, the stabilizing effects of viscosity and thermal conductivity are overcome by the destabilizing buoyancy, and an overturning instability ensues as thermal convection: hotter part of fluid is lighter and tends to rise as colder part tends to sink according to the action of the gravity force [12]. The phenomenon of thermal convection itself had been recognized by Rumford [60] and Thomson [70]. However, the first quantitative experiment on thermal instability and the recognition of the role of viscosity in the phenomenon are due to Bénard [1], so the convection in a horizontal layer of a fluid heated from below is called Bénard convection. For many years, the question for understanding of convective flows has motivated numerous theoretical, numerical, and experimental studies [22,12].

Thermal convection in viscoelastic fluids is also a subject of considerable interest in contemporary fluid flow and heat transfer researches. The first work which deals directly with thermal instability of a viscoelastic fluid appears to be that of Herbert who studied plane Couette flow heated from below [21]. He found that a finite elastic stress in the undisturbed state is necessary for elasticity to affect stability. Since Herbert's pioneering work, many physicists have continued to develop the linear theory and nonlinear numerical method in the studies of thermal instability in viscoelastic fluids, see [11,14,47,40,59,61,65,66] and the references cited therein. Moreover, it has also been widely investigated how thermal convection in viscoelastic fluids evolves under the effects of other physical factors, such that rotation [13,42,67], magnetic fields [2,3,54], the porous media [77,62] and so on.

As Rosenblat pointed out in [59], the nature of (linear) convective solution depends strongly on the particular constitutive relation used to characterize the viscoelasticity. For certain models and certain parameter ranges the convection is supercritical and stable, while for other models and parameter ranges it can be subcritical and unstable. In other words, the influence of elasticity on the thermal convection is closely related to the choice of model describing the motion of a viscoelastic fluid.

Recently Jiang–Liu mathematically prove the phenomenon of inhibition of thermal instability by elasticity by the following (nonlinear) Boussinesq approximation equations of viscoelastic fluids [38]:

We shall explain the mathematical notations in (3.1). The unknowns $v: = v(x,t)$, $\Theta: = \Theta(x,t)$, $U: = U(x,t)$ and $p: = p(x,t)$ represent the velocity, temperature, deformation tensor and pressure of an incompressible viscoelastic fluid resp.. The parameters $\rho$, $\alpha$, $\kappa$ and $g>0$ denote the density constant at some properly chosen temperature parameter $\Theta_{ \mathrm{b}}$, the coefficient of volume expansion, the elasticity coefficient of the fluid and the gravitational constant resp.. $k: = \sigma/\rho c_{ \mathrm{v}}$ and $\nu: = \mu/\rho$ represent the coefficient of the thermometric conductivity and the kinematic viscosity resp., where $c_{ \mathrm{v}}$ is the specific heat at constant volume, $\sigma$ the coefficient of heat conductivity and $\mu$ the coefficient of shear viscosity. $e_{3} = (0,0,1)^{ \mathrm{T}}$ stands for the vertical unit vector, $g\rho\alpha(\Theta-\Theta_{ \mathrm{b}})e_{3}$ for the buoyancy (caused by expanding with heat and contracting with cold) and $-\rho g e_{3}$ for the gravitational force.

The rest state of the Boussinesq approximation equations (3.1) can be given by $r_{ \mathrm{V}}: = (0,\bar{\Theta},I)$ with an associated pressure profile $\bar{p}$, where the temperature profile $\bar{\Theta}$ and $\bar{p}$ depend on $x_{3}$ only, and satisfy the equilibrium (or rest) state:

For the simplicity, we consider that

where $\varpi>0$ is a constant of adverse temperature gradient. It should be noted that in this section we define that

where $\mathbb{T}$ is defined by (2.2).

Denoting the perturbation around the equilibrium state by

then, $(v,\theta,V,\beta)$ satisfies the following perturbation system of equations

We shall pose the following initial-boundary value conditions for the well-posedness of (3.3):

We call the initial-boundary value problem (3.3)–(3.5) the viscoelastic Rayleigh–Bénard problem.

Now we set the unknowns in Lagrangian coordinates by

where $\zeta$ is defined by (1.2) with $\Omega_h$ being in place of $\Omega$. We also assume that $\det(\nabla \eta^0+I)^{- \mathrm{T}} = 1$ and (1.3) are satisfied in $\Omega_h$. Thus, the evolution equations for $(u,\vartheta,q)$ in Lagrangian coordinates read as follows

with initial-boundary value conditions

From now on, we still define that $\eta: = \zeta-y$.

Let $Q = {h^{2}\kappa}/{\rho\nu^{2}}$, $R^{2} = {g\varpi\alpha h^{4}}/{\nu k}$ and $P_{\vartheta} = {\nu}/{k}$, where $Q$, $R^{2}$ and $P_{\vartheta}$ are called the elasticity, Rayleigh and Pandtl numbers, resp.. Now we use dimensionless variables

to rewrite (3.6) as the following non-dimensional form defined in $\Omega$:

where $\mathbb{D}\eta: = \nabla\eta+\nabla\eta^{ \mathrm{T}}$. The corresponding initial-boundary value conditions read as follows

We call the initial-boundary value problem (3.7)–(3.8) the transformed VRB problem. Before stating the main results for the transformed VRB problem (2.7), we shall introduce some notations of Sobolev spaces in this section: for $k\geqslant 1$,

We have the following stability result for the transformed VRB problem:

Theorem 3.1. [38,Theorem 2.1] Under the condition

then there is a sufficiently small $\delta>0$, such that for any $(\eta^{0},u^{0},\vartheta^{0})\in H_{0}^{3}(\Omega)\times H_{0}^{2}(\Omega)\times H_{0}^{2}(\Omega)$ satisfying $\det(\nabla \eta^0+I) = 1$, $\mathrm{div}_{\mathcal{A}^0}u^0 = 0$ and

there exists a unique strong solution $(\eta,u,\vartheta,q)\in C^{0}([0,\infty),H_{0}^{3}(\Omega)\times H^{2}_0(\Omega)\times H^{2}(\Omega)\times \underline{H}^1(\Omega))$ to the transformed VRB problem (3.7)–(3.8). Moreover, $(\eta,u,\vartheta,q)$ enjoys the following exponential stability estimate:

The above positive constants $\delta$, $c$ and $\omega$ depend on the domain $\Omega$ and other known physical parameters in the transformed VRB problem.

The viscoelastic Rayleigh–Bénard problem (3.3)–(3.5) in the absence of deformation tensor reduces to the classical RB (i.e., Rayleigh–Bénard) problem

We mention that there exists a threshold $R_0$ defined by

such that if the convection condition $R>R_0$ is satisfied, the RB problem is unstable; if $R<R_0$, the RB problem is stable. Hence the stability result in Theorem 3.1 presents that the elasticity can inhibit the thermal instability for sufficiently large $\kappa$.

Next we briefly explain why the stability condition is given by the form (3.9). Let us first recall the basic energy identity:

By the idea of the inhibition of instability by the elasticity, (3.7)$_1$ and (3.7)$_3$, we rewrite the above identity as follows:

where

To control the first integral on the right hand of (3.10), we shall derive from (3.7)$_1$ and (3.7)$_2$ that for any multiindex $\alpha = (\alpha_1,\alpha_2)$ satisfying $0\leqslant |\alpha | \leqslant 1$,

Here and in what follows, $\partial_{ \mathrm{h}}^\alpha : = \partial_1^{\alpha_1} \partial_2^{\alpha_2}$. Hence, we further derive from the above two identities that

where

It is easy to see that, for sufficiently large $Q$,

Moreover precisely, they are equivalent to, for sufficiently large $Q$,

resp.. Based on the above idea, we can easily find out the stability condition (3.9), which is relatively more complicated than the stability condition (2.8) in the transformed VRT problem (2.7). Moreover, we can establish Theorem 3.1 under the stability condition by an energy method.

In addition, we also find that the thermal instability can occur if $\kappa$ is sufficiently small. More precisely, we have the following instability result.

Theorem 3.2. [37,Theorem 2.2] We define that

Let $R$ satisfy the convection condition $R>R_{0}$. If $Q$ and $R$ satisfy the condition

where we have defined that $H: = R/\sqrt{P_\vartheta}-2(R-R_0)\xi$, the zero solution to the transformed VRB problem (3.7)–(3.8) is unstable in the Hadamard sense, that is, there are positive constants $\Lambda^{*}$, $m_{0}$, $\epsilon$, $\delta_{0}$ and

such that, for any $\delta\in(0,\delta_{0})$ and the initial data

there is a unique strong solution $(\eta,(u,\vartheta) )\in C^{0}([0,T],H_{0}^{3}\times H_{0}^{2},\underline{H}^{1}(\Omega) )$ with initial data $(\eta^{0},u^{0},\vartheta^{0})$ to the transformed VRB problem. However the solution satisfies (2.13) for some escape time $T^{\delta}: = \frac{1}{\Lambda^{*}}\ln\frac{2\epsilon}{m_{0}\delta}\in I_{T}$. Moreover, the initial data $(\eta_{0},u_{0} )$ satisfies $\det(\nabla \eta^0+I) = 1$ and $\mathrm{div}_{\mathcal{A}^{0}}u_{0} = 0$ in $\Omega$, where $\mathcal{A}^{0}: = (\nabla \eta^0+I)^{- \mathrm{T}}$.

Remark 8. Similarly to (6), if $\delta$ is sufficiently small, the solution obtained in Theorems 2.1 and 2.2 automatically satisfies the properties:

Thus we can also recover the exponential stability of the VRB problem (3.3)–(3.5) in Eulerian coordinates from Theorems 2.1 and 2.2 by an inverse transformation.

Similarly to Theorem 2.2 we can use a bootstrap instability method in [16] to establish Theorem 3.2. However, the construction of unstable linear solutions is relatively more complicated. In fact, following the argument of linear instability of the linearized VRT problem in Remark 7, we will face the following variational problem

where $\mathcal{A}$ is defined by (3.11). In particular, we see that the stabilizing term $Q\|\nabla u\|_{{L^2(\Omega)}}^2$ to inhibit the thermal instability term $- 2R\int_\Omega u_3\vartheta\mathrm{d}y$ depends on $\lambda$. This dependence relation results in the complexity of the above variational problem, and the difficulty of looking for the threshold value of $Q$ for the instability. However we can find out an instability condition (3.12) by careful analysis.

4.   Existence of large solutions of some class

The two stability results in the previous sections present that the elasticity has the stabilizing effect in the motion of viscoelastic fluids under the small perturbations. In this section we will consider that this stabilizing effect also be observed under the large perturbations.

From now on, we define $\mathbb{T}: = \mathbb{T}_1\times \mathbb{T}_2\times \mathbb{T}_3$, where $\mathbb{T}_i = 2\pi L_i(\mathbb{R}/\mathbb{Z})$, and $2\pi L_i$ ($i = 1$, $2$, $3$) are the periodicity lengths. Let $( u ,U, q)$ satisfy the system (1.1) with initial data $(v^0,U^0)$ in $\mathbb{T}$ and $\zeta$ be defined by (1.2) with $\mathbb{T}$ being in place of $\Omega$. We also assume that (1.3) is satisfied in $\mathbb{T}$. Let $\eta: = \zeta-y$ and

Then the evolution equations for $(\eta,u,q)$ in Lagrangian coordinates read as follows.

with initial data

Now we further define some simplified function spaces:

Here and in what follows, $(w)_{\mathbb{T}} = \int_{\mathbb{T}}w \mathrm{d}y/8L_1L_2L_3\pi^3$ for $w\in L^1(\mathbb{T})$.

Next we introduce the first result, which is concerned with the existence of strong solutions to the initial value problem (4.1)–(4.2) in some classes of large initial data:

Theorem 4.1. [33,Theorem 1.2] There are constants ${c}_1\geqslant 1$ and ${c}_2\in (0,1]$, such that for any $(\eta^0,u^0)\in H^{3}_{*,1}(\mathbb{T} )\times H^{2} (\mathbb{T} )$ and $\kappa$, satisfying $\mathrm{div}_{\mathcal{A}^0}u^0 = 0$ and

where $\mathcal{A}^{0}: = (\nabla \eta^0+I)^{- \mathrm{T}}$ and $I_0^{ \mathrm{h}}( u^0,\eta^0): = \|(u^0,(1+\sqrt{\kappa}) \nabla \eta^0) \|_{H^2(\mathbb{T} )}^2$, the initial value problem (4.1)–(4.2) admits a unique global solution $(\eta, u,q)\in C^0([0,\infty), H^3_{*,1}(\mathbb{T} ))\times C^0([0,\infty)$, $H^2(\mathbb{T} ))\times C^0([0,\infty), \underline{H}^2(\mathbb{T} ) )$. Moreover, the solution $(\eta,u)$ enjoys the stability estimate

Here and in what follows the positive constants $c$, $c_1$–$c_3$ depend on $\rho$ and $\mu$, and $c$ varies from line to line.

Theorem 4.1 exhibits the global existence of strong solutions to the initial value problem of (1.1) defined in a spatially periodic domain, when the initial deformation (i.e., $U-I$ at $t = 0$) and the initial velocity are small for given parameters. However the initial velocity can be large if the elasticity coefficient is large. This means that the strong elasticity can prevents the development of singularities even when the initial velocity is large, thus playing a similar role to viscosity in preventing the formation of singularities in pure viscous flows.

We briefly sketch the proof idea of Theorem 4.1. Recalling (1.10), we easily conclude that the equation (4.1) may be approximated by the following linear system of equations for sufficiently large $\kappa$:

Since the linear system of equations has global solutions with large initial data, we could expect that the initial value problem (4.1)–(4.2) may also admit a global large solution for sufficiently large $\kappa$ (with fixed $I_0^{ \mathrm{h}}(u^0,\eta^0)$) as stated in Theorem 4.1.

In order to obtain Theorem 4.1, the key step is to derive the a priori estimate (4.3) under sufficiently large $\kappa$. Fortunately, by careful energy estimates, we find that the estimate (1.10) is still inherited by the higher order spatial derivatives of $(\eta,u)$ under sufficiently large $\kappa$. Namely, we can conclude that there are constants $K$ (depending possibly on $I_0^{ \mathrm{h}}(u^0,\eta^0)$ but not on $T$) and $\delta$, such that

provided that

and

Based on the above fact and the existence of a unique local solution, we can immediately obtain Theorem 4.1.

Next we further introduce the second result, which is concerned with the properties of the solution $(\eta,u)$ to the initial value problem (4.1)–(4.2) in Theorem 4.1. In particular, we will see that the motion of the viscoelastic fluid can be approximated by a linear pressureless motion in Lagrangian coordinates for sufficiently large $\kappa$, even when the initial velocity is large.

Theorem 4.2. [33,Theorem 1.2] Let $(\eta,u)$ be the solution established in Theorem 4.1, then we have

(1) Exponential stability of $(\eta,u)$: for any $t\geqslant 0$,

where $\bar{u}: = u- (u)_{\mathbb{T}}$ and $\bar{u}^0: = u^0-(u^0)_{\mathbb{T}}$.

(2) Large-time behavior of $\eta$:

where $\bar{\eta}: = \eta-\eta(y_1,y_2,0 ,t)$ and $\varpi = (\eta^0)_{\mathbb{T}}$.

(3) Stability of $(\eta,u)$ around $( \eta^{ \mathrm{l}}, u^{ \mathrm{l}})$:

Here $(\eta^{ \mathrm{d}}, u^{ \mathrm{d}} ): = (\eta-\eta^{ \mathrm{l}} , u-u^{ \mathrm{l}} )$ and $(\eta^{ \mathrm{l}}, u^{ \mathrm{l}} )\in C^0([0,\infty), H^{3} (\mathbb{T} ))\times C^0([0,\infty)$, $H^{2} (\mathbb{T} ))$ is the unique strong solution of the linear system (4.4) with $q^{ \mathrm{l}} = 0$ and initial data:

where $(\eta^{ \mathrm{r}}, u^{ \mathrm{r}})\in \underline{H}^3 (\mathbb{T})\times \underline{H}^2 (\mathbb{T})$ has the following properties:

(a) $\mathrm{div}(u^0+u^{ \mathrm{r}}) = 0$ and $\|u^{ \mathrm{r}} \|_{H^2(\mathbb{T} )}\leqslant c_0 \|\nabla \eta^0\|_{H^2(\mathbb{T})}\|u^0\|_{H^2(\mathbb{T} )}$.

(b) $\mathrm{div}(\eta^0+\eta^{ \mathrm{r}}) = 0$ and $\|\eta^{ \mathrm{r}} \|_{H^3(\mathbb{T} )}\leqslant c_0 \|\nabla \eta^0\|_{H^2(\mathbb{T} )}^2$,

where the constants $c_0$ is independent of any paraments.

We end this section by listing some remarks for Theorem 4.2.

Remark 9. We explain on why the initial data $(\eta^0,u^0)$ has to be modified as in (4.8).

(1) Since the initial data for $u^{ \mathrm{l}}$ has to satisfy the divergence-free condition, i.e., $\mathrm{div}(u^{ \mathrm{l}}|_{t = 0}) = 0$, one thus has to adjust the initial data $u^0$ as in (4.8).

(2) The initial data $\eta^0$ for $\eta$ can be directly used as an initial data for $\eta^1$. In this case one can see that $\mathrm{div}\eta^{ \mathrm{l}} = \mathrm{div}\eta^0$. Consequently, the decay-in-time of $\nabla \eta^{ \mathrm{d}}$ by (4.5) can not be expected, unless $\mathrm{div}\eta^0 = 0$. Hence, we have to modify $\eta^0$ as in (4.8), so that the new initial data"$\eta^0+\eta^{ \mathrm{r}}$" also satisfies the divergence-free condition.

Remark 10. Let us try to give a physical meaning hidden in (4.6). Consider the viscoelastic fluid in a periodic cell $K: = (0,2\pi L_1)\times(0,2\pi L_2)\times(0,2\pi L_3)$, and think that the fluid is made up of infinite fluid segments that are parallel to $x_3$-axis. Now, we consider a straight line segment denoted by $l^0$ in $K$ at $t = 0$, and any given (fluid) particle $y$ in $l^0$. The two particles at the upper and lower endpoints are denoted by $y^1$ and $y^2$, and note that $y_{ \mathrm{h}} = y_{ \mathrm{h}}^1 = y_{ \mathrm{h}}^2$. We disturb the rest state by a perturbation $(\eta^0,u^0)$ at $t = 0$. Then, the line segment $l^0$ will be bent and move to a new location at time $t$ that may be curved line segment and denoted by $l^t$. Since the motion of the fluid is spatially periodic, the segment $y^1y^2$ is thus parallel to $x_3$-axis and can be represented by

At time $t$, the new location of a particle $y$ on $l^0$ is given by $\eta(y,t)+y$. Thus we see that $|\bar{\eta} (y,t)|$ has a geometric meaning, namely, it represents the distance between the two points $\eta(y,t)+y$ on $l^t$ and $(\eta_{ \mathrm{h}}(y_{ \mathrm{h}},0,t)+y_{ \mathrm{h}}, \eta_3(y_{ \mathrm{h}},0,t)+y_3)$ on $l^{ \mathrm{n}}$.

By (4.6) and the interpolation inequality, we see that

In particular, for the case $\eta^0 = 0$,

from which and the geometric meaning of $|\bar{\eta} (y,t)|$ we immediately see that the curve $l^t$ oscillates around $l^{ \mathrm{n}}$, and the amplitude tends to zero, as $\kappa$ or $t$ is extremely large.

Remark 11. Similarly to Remark 10, we can also give a physical meaning hidden in (4.7). Namely, consider a straight line segment $l^0$ consisted of fluid particles in the rest state and perturb the segment by a velocity, thus it will be bent. Then, as $\kappa$ or $t\to\infty$, the perturbed segment will turn into a straight line segment that is parallet to $l^0$ and has the same length as $l^0$.

Remark 12. The corresponding versions of Theorems 4.1 and 4.2 can be found in Theorems 1.4 and 1.5 in [33].

5.   Stabilizing effect in elasticity fluids

In previous three sections, we have seen the stabilizing effect of elasticity on the motion of viscoelastic fluids. In fact, this stabilizing effect has been also mathematically verified in the elasticity fluids. Next we further briefly introduce the relevant results in elastic fluids.

The system without viscosity reduces to the following system, which describes the motion of elasticity fluids:

In the absence of $U$ in (5.1), we obtain the famous Euler equations for the motion of an inviscid, incompressible fluid:

The global regularity of solutions to the two-dimensional Euler equations has been known for a long time. It is also known that the gradient of the vorticity given by $\mathrm{curl}v$ and its higher-order Sobolev norms may grow unboundedly. The best known upper bound on the growth of the gradient of vorticity is double exponential in time. Kiselev–$\mathrm{\check{S}}$verák had constructed an example of initial data in the disk such that the corresponding solution to the 2D Euler equations exhibits double exponential growth in the gradient of vorticity for all times [41]. Later Zlato$\mathrm{\check{s}}$ further showed that in the 2D periodic domain the vorticity gradient can grow at least exponentially as $t\to \infty$ [80]. At present the analogous question on the whole space is still open.

However Lei proved the existence of global solutions with small initial data to the 2D Cauchy problem of (5.1) in Lagrangian coordinates, where the highest-order energy solutions at most algebraically grows in time [44]. Recently the uniform bound of the highest-order energy of global solutions to the 2D case was further proved by Cai [6], which exhibits that the elasticity can inhibit the growth of solutions. We mention that Wang also gave the existence of global solutions to the 2D Cauchy problem of (5.1) by using space-time resonance method and a normal form transformation [72], and the relevant result of the 3D system of (5.1) can be found in [64].

In addition, the stabilizing effect of elasticity on the local-in-time motion of elasticity fluids can be found in the free-boundary case, interested readers can refer to [8] for the vortex sheet problem and [46] for the RT problem.

6.   Conclusions

In this review, we have introduced the mathematical results concerning the stabilizing effect of elasticity on the motion of viscoelastic/elastic fluids. However these results were verified by the Oldroyd model, which includes a viscous stress component and a stress component for a neo-Hookean solid. We except that similar stabilizing results can be extended to the other Oldroyd-B models in future. In addition, there are still many relevant interesting stabilizing problems, which should be further investigated, for examples,

(1) Theorem 2.1 provides the sharp stability criterion for the transformed VRT problem. What about the critical case $\mathrm{C}_{ \mathrm{r}} = 1$? Under such case, we physically guess the transformed VRT problem is stable for some initial data, and unstable for other initial data by the minimal potential energy principle.

(2) Theorem 3.1 provides the stability condition of $Q$ for the viscoelastic Rayleigh–Bénard problem. It still is extremely difficult to find out the threshold value of $Q$.

(3) Theorem 4.1 proves the existence of large solutions of some class in the spatially periodic domain case. It is not clear that whether a similar result can be found in a general bounded domain.

Acknowledgments

The author is grateful to Dr. Binqiang Xie in South China Normal University for communication in the vortex sheet problem of elasticity fluids.