On multiple solutions to the steady flow of incompressible fluids subject to do-nothing or constant traction boundary conditions on artificial boundaries

TL;DR

This paper investigates the implications of using do-nothing boundary conditions in steady incompressible fluid flow simulations, revealing scenarios with multiple solutions or solution blow-up, highlighting issues with well-posedness.

Contribution

It provides analytical and numerical examples demonstrating multiple solutions and blow-up phenomena caused by do-nothing boundary conditions in incompressible flow problems.

Findings

Unsteady solutions can blow up or converge depending on initial data.

Multiple steady solutions can exist under do-nothing boundary conditions.

Numerical simulations confirm analytical findings in complex geometries.

Abstract

The boundary conditions prescribing the constant traction or the so-called do-nothing conditions are frequently taken on artificial boundaries in the numerical simulations of steady flow of incompressible fluids, despite the fact that they do not facilitate a well-posed problem, not allowing to establish the standard energy estimates. In a pursuit to understand better the possible consequences of using these conditions, we present a particular set of examples of flow problems, where we find none or two analytical or numerical solutions. Namely, we consider problems where the flow connects two such artificial boundaries. In the simple case of the isotropic radial flows where both steady and unsteady analytical solutions are derived easily, we demonstrate that while for some (large) boundary data all unsteady solutions blow up in finite time, for other data (including small or trivial) the unsteady flows either converge asymptotically to one of two steady solutions, or blow up in finite time, depending on the initial state. We then document the very same behavior of the numerical solutions for planar flow in a diverging channel. Finally, we provide an illustrative example of two steady numerical solutions to the flow in a three-dimensional bifurcating tube, where the inflow velocity is prescribed at the inlet, while the two outlets are treated by the do-nothing boundary condition.

Equations

$$\left.\begin{array}[]{rcl}\operatorname{div}(\boldsymbol{u}\otimes\boldsymbol{u})-\operatorname{div}\mathbf{T}&=&\boldsymbol{f}\\ \operatorname{div}\boldsymbol{u}&=&0\end{array}\quad\right\}\quad\textrm{ in $\Omega$, where}\quad\begin{array}[]{rcl}\mathbf{T}&=&-p\mathbf{I}+\mathbf{S},\\ \mathbf{S}&=&\nu\,(\nabla\boldsymbol{u}+(\nabla\boldsymbol{u})^{T}).\end{array}$$

$$(\mathrm{div}(\mathbf{u}\otimes \mathbf{u}),\mathbf{w}{)}_{\Omega }-(\mathrm{div}\mathbf{T},\mathbf{w}{)}_{\Omega }$$

$$(\mathrm{div}\mathbf{u},q{)}_{\Omega }$$

$$(\mathrm{div}\mathbf{T},\mathbf{w}{)}_{\Omega }=(\mathbf{T}\mathbf{n},\mathbf{w}{)}_{\partial \Omega }-(\mathbf{T},\nabla \mathbf{w}{)}_{\Omega },$$

$$\mathbf{u}=0\text{on Γ0⊂∂Ω},$$

$$\begin{array}[]{rcl}-\mathbf{T}\boldsymbol{n}=p\boldsymbol{n}-\mathbf{S}\boldsymbol{n}=&&\\ p\boldsymbol{n}-\nu\,(\nabla\boldsymbol{u}+(\nabla\boldsymbol{u})^{T})\boldsymbol{n}&=&\boldsymbol{b}\qquad\textrm{on $\Gamma_{\boldsymbol{b}}\subset\partial\Omega$},\end{array}$$

$$(\mathrm{div}(\mathbf{u}\otimes \mathbf{u}),\mathbf{w}{)}_{\Omega }+(\nu (\nabla \mathbf{u}+(\nabla \mathbf{u}{)}^T),\nabla \mathbf{w}{)}_{\Omega }-(p,\mathrm{div}\mathbf{w}{)}_{\Omega }-(q,\mathrm{div}\mathbf{u}{)}_{\Omega }=(\mathbf{f},\mathbf{w}{)}_{\Omega }-(\mathbf{b},\mathbf{w}{)}_{{\Gamma }_{\mathbf{b}}}.$$

$$\mathrm{div}(\nu (\nabla \mathbf{u}+(\nabla \mathbf{u}{)}^T))=\nu \mathrm{div}(\nabla \mathbf{u})+\nu \nabla (\mathrm{div}\mathbf{u})=\nu \mathrm{div}(\nabla \mathbf{u})=\nu \Delta \mathbf{u}.$$

$$\mathrm{div}(\mathbf{u}\otimes \mathbf{u})-\nu \Delta \mathbf{u}+\nabla p=\mathbf{f}\text{in Ω}.$$

$$\hat{\mathbf{T}}=-p\mathbf{I}+\hat{\mathbf{S}},\hat{\mathbf{S}}=\nu \nabla \mathbf{u},$$

$$\begin{array}[]{rcl}-\hat{\mathbf{T}}\boldsymbol{n}=p\boldsymbol{n}-\hat{\mathbf{S}}\boldsymbol{n}=&&\\ p\boldsymbol{n}-\nu\,(\nabla\boldsymbol{u})\boldsymbol{n}&=&\hat{\boldsymbol{b}}\qquad\textrm{on $\Gamma_{\hat{\boldsymbol{b}}}\subset\partial\Omega$}.\end{array}$$

$$(\mathrm{div}(\mathbf{u}\otimes \mathbf{u}),\mathbf{w}{)}_{\Omega }+(\nu \nabla \mathbf{u},\nabla \mathbf{w}{)}_{\Omega }-(p,\mathrm{div}\mathbf{w}{)}_{\Omega }-(q,\mathrm{div}\mathbf{u}{)}_{\Omega }=(\mathbf{f},\mathbf{w}{)}_{\Omega }-(\hat{\mathbf{b}},\mathbf{w}{)}_{{\Gamma }_{\hat{\mathbf{b}}}}.$$

$$(\mathrm{div}(\mathbf{u}\otimes \mathbf{u}),\mathbf{u}{)}_{\Omega }+{\left(\nu (\nabla \mathbf{u}+(\nabla \mathbf{u}{)}^T),\frac{\nabla \mathbf{u}+(\nabla \mathbf{u}{)}^T}{2}\right)}_{\Omega }$$

$$\text{or}(\mathrm{div}(\mathbf{u}\otimes \mathbf{u}),\mathbf{u}{)}_{\Omega }+(\nu \nabla \mathbf{u},\nabla \mathbf{u}{)}_{\Omega }$$

$$(\mathrm{div}(\mathbf{u}\otimes \mathbf{u}),\mathbf{u}{)}_{\Omega }+\sigma \Vert \tilde{\mathbf{S}}{\Vert }_{2,\Omega }^2\le \Vert \mathbf{f}{\Vert }_{2,\Omega }\Vert \mathbf{u}{\Vert }_{2,\Omega }+\Vert \tilde{\mathbf{b}}{\Vert }_{2,{\Gamma }_{\tilde{\mathbf{b}}}}\Vert \mathbf{u}{\Vert }_{2,{\Gamma }_{\tilde{\mathbf{b}}}},$$

$$(\mathrm{div}(\mathbf{u}\otimes \mathbf{u}),\mathbf{u}{)}_{\Omega }=\frac{1}{2}(\mathrm{div}\mathbf{u},|\mathbf{u}{|}^2{)}_{\Omega }+\frac{1}{2}(\mathbf{u}\cdot \mathbf{n},|\mathbf{u}{|}^2{)}_{\partial \Omega }=\frac{1}{2}(\mathbf{u}\cdot \mathbf{n},|\mathbf{u}{|}^2{)}_{{\Gamma }_{\tilde{\mathbf{b}}}},$$

$${\int }_{{\Gamma }_{\tilde{\mathbf{b}}}}\tilde{\mathbf{b}}(\mathbf{u})\cdot \mathbf{u}\mathrm{d}S\ge -\frac{1}{2}{\int }_{{\Gamma }_{\tilde{\mathbf{b}}}}(\mathbf{u}\cdot \mathbf{n})|\mathbf{u}{|}^2\mathrm{d}S-\text{ lower order terms}$$

$$\tilde{\mathbf{b}}(\mathbf{u})=P_0\mathbf{n}-\frac{1}{2}(\mathbf{u}\cdot \mathbf{n})\mathbf{u}\text{or}\tilde{\mathbf{b}}(\mathbf{u})=\left(P_0-\frac{1}{2}|\mathbf{u}{|}^2\right)\mathbf{n}$$

$$\hat{\mathbf{b}}(\mathbf{x},\mathbf{u})=\mathbf{h}(\mathbf{x})+\alpha (\mathbf{u}\cdot \mathbf{n}{)}^{-}\mathbf{u}\text{for }\alpha >\frac{1}{2},$$

$$\hat{\mathbf{b}}=P_0\mathbf{n},\text{i.e.}(p-P_0)\mathbf{n}-\nu (\nabla \mathbf{u})\mathbf{n}=0,$$

$$\tilde{\mathbf{b}}=P_i\mathbf{n}\text{on }{\Gamma }_i,i=1,2,$$

$$\mathbf{u}=U(r)\mathbf{x},\text{with }r=|\mathbf{x}|,\mathbf{x}=(x,y).$$

$$0=\mathrm{div}\mathbf{u}=2U+r{U}^{\prime }⟹U(r)=r^{-2}Q,$$

$$\mathrm{div}(\mathbf{u}\otimes \mathbf{u})=(3U^2+2rU{U}^{\prime })\mathbf{x}\stackrel{\text{eqref}rad:div}{=}-r^{-4}{Q}^2\mathbf{x}.$$

$$\nabla\boldsymbol{u}=(\nabla\boldsymbol{u})^{T}=\frac{1}{r}\left(\begin{array}[]{cc}rU+x^{2}U^{\prime}&xyU^{\prime}\\ xyU^{\prime}&rU+y^{2}U^{\prime}\end{array}\right)\overset{\eqref{rad:div}}{=}\frac{Q}{r^{4}}\left(\begin{array}[]{cc}r^{2}-2x^{2}&-2xy\\ -2xy&r^{2}-2y^{2}\end{array}\right),$$

$$\mathrm{div}(\nu \nabla \mathbf{u})=\mathrm{div}(\nu (\nabla \mathbf{u}{)}^T)=\nu (3r^{-1}{U}^{\prime }+U^{\prime \prime })\mathbf{x}\stackrel{\text{eqref}rad:div}{=}0.$$

$$p=P(r),\text{so that}\nabla p=r^{-1}{P}^{\prime }\mathbf{x},$$

$$0=(3U^2+2rU{U}^{\prime })-2\nu (3r^{-1}{U}^{\prime }+U^{\prime \prime })+r^{-1}{P}^{\prime }\stackrel{\text{eqref}rad:div}{=}-r^{-4}{Q}^2+r^{-1}{P}^{\prime },$$

$$P^{\prime }=r^{-3}{Q}^2⟹P=P_{\infty }-\frac{1}{2}{r}^{-2}{Q}^2,\text{with some }{P}_{\infty }\in \mathbb{R}.$$

$$\begin{array}{rl}\Omega & =\{R_1^2<x^2+y^2<R_2^2,x>0,0<y/x<\tan (\alpha )\},\\ {\Gamma }_i& =\partial \Omega \cap \{x^2+y^2=R_i^2\},i=1,2,\\ {\Gamma }_{\mathrm{n}}& =\partial \Omega \cap (\{x=y/\tan (\alpha )\}\cup \{y=0\}).\end{array}$$

$$\mathbf{u}\cdot \mathbf{n}=0\text{and}(-\mathbf{T}\mathbf{n}{)}_{\tau }=0\text{on }{\Gamma }_{\mathrm{n}}.$$

Full text

1 Introduction

The most of studies in fluid dynamics or computational fluid dynamics deal with flows in domains smaller than the whole space $\mathbb{R}^{2}$ or $\mathbb{R}^{3}$. The governing equations of the flow then consist of a system of partial differential equations inside the domain $\Omega$ and additional conditions given on the boundary $\partial\Omega$, only together can they define a well-posed problem. On the parts of the boundary that are related to a natural material interface, the physics of the interface is employed to derive the boundary conditions, similarly as when deriving the partial differential equations in the bulk based on the rheology of the fluid. The no-slip assumption considered on the interface of the fluid with a solid wall and the corresponding Dirichlet boundary conditions can serve as an elementary example, one among others. In many practical problems however, some parts of the boundary are not related to any such interface but represent a truncation of the flow extending beyond the considered domain; they are artificial boundaries. The boundary conditions on artificial boundaries cannot be derived merely from the underlying physics, but are justified by a combination of arguments related to the physics, the modelling goals, the well-posedness analysis and the numerical analysis.

This paper is concerned with two types of boundary conditions that are frequently used in numerical simulations of flows of incompressible fluids, in particular Navier–Stokes fluid, despite the fact that they do not facilitate the well-posedness of the problem. Namely, as introduced in a greater detail below, the boundary conditions prescribing a constant traction, or the ones that we later call the full-gradient-traction conditions and cover the so-called do-nothing conditions, are considered. It is well known that the steady Navier–Stokes equations subject to these boundary conditions do not allow for standard energy estimates. Consequently, it has been impossible so far to establish the existence theory for the steady problem except for small data (and similarly, to establish the existence of unsteady solution except for small data or small time interval). While the non-uniqueness of steady solutions seems to be expected intuitively, no concrete example of multiple solutions has been given in the literature so far, to the best of our knowledge. Indeed, as Galdi comments in [1], p. 180, (after proving that for small data there is a unique small solution)

“…the question of whether a given solution is unique in the class of all possible weak solutions corresponding to the same data …is, to date, open, in the case of the do-nothing conditions.”

and similarly does Rannacher in [2], p. 294,

“Unfortunately, in the case of open outlets such an a priori bound is not known. This is reflected by the fact that not even the global uniqueness of the zero-solution has been proven yet …However, this possible non-uniqueness could not be confirmed by numerical experiments.”

The very last sentence of the above quotation points out an important aspect: that, in fact, the addressed boundary conditions are used in numerical simulations by many researchers on a regular basis, since in many problem geometries and flow regimes they are experienced to deliver a unique solution.

This paper is motivated by the above stated incongruity between the hitherto achievements in the theoretical analysis of Navier–Stokes equations, and the well-established techniques and needs in practical numerical simulations. Our aim is to contribute to a better understanding of the corresponding solutions, or even to provoke further study of the theoretical aspects of Navier–Stokes equations in this respect. The scope of the paper is to present a set of examples of multiple solutions to the steady flow subject to the above addressed boundary conditions. Aside from the mere confirmation of the non-uniqueness of solutions we intend to demonstrate a particular simple mechanism behind, and to pursue the behavior of the steady solutions, of the corresponding unsteady flows and of their numerical approximation. We observe multiple solutions for small boundary data, including the case of trivial data where both the trivial and non-trivial solutions can be found. On the other hand, for some instances of large boundary data, the considerations in a simplified setting and our numerical simulations indicate the possibility that no steady solution would exist, this fact being related to the very same mechanism.

1.1 Weak solutions, traction, full-gradient-traction, and do-nothing

The steady Navier–Stokes equations can be written as

[TABLE]

Here $p$, $\boldsymbol{u}$ denote the kinematic pressure and the velocity, $\mathbf{T}$ and $\mathbf{S}$ stand for Cauchy stress tensor and its viscous part, the constant $\nu>0$ being the kinematic viscosity, and $\boldsymbol{f}$ represents the density of outer forces (later on, we set $\boldsymbol{f}\equiv\boldsymbol{0}$ for simplicity). The domain $\Omega\subset\mathbb{R}^{2}$ or $\Omega\subset\mathbb{R}^{3}$ is considered bounded and with Lipschitz boundary $\partial\Omega\in\mathcal{C}^{0,1}$.

The weak formulation can be formally derived by multiplying the system (1) by suitable test functions $\boldsymbol{w}$, $q$ and then integrating over $\Omega$. One arrives at

[TABLE]

to hold for all $\boldsymbol{w}$, $q$. The suitable function spaces for $\boldsymbol{u}$, $\boldsymbol{w}$ and $p$, $q$ are to be specified depending on the boundary conditions considered. By the brackets we shortly denote the integrals of the scalar product $(\boldsymbol{f},\boldsymbol{g})_{\Omega}=\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{g}\,\mathrm{d}\boldsymbol{x}$, and analogously for scalars or tensors and for the integration on $\partial\Omega$. The divergence theorem allows us to write

[TABLE]

where $\boldsymbol{n}$ denotes the unit outer normal vector. From this it appears that the weak formulation offers two quantities to be naturally considered on the boundary: the velocity $\boldsymbol{u}$ (to be treated by an essential, Dirichlet, boundary condition) and the surface traction $-\mathbf{T}\boldsymbol{n}$ (to be imposed by a natural, Neumann, boundary condition). Imposing, for instance, the no-slip conditions on a part of the boundary

[TABLE]

and prescribing the traction $\boldsymbol{b}$ on the remaining part (the artificial boundary)

[TABLE]

the weak solution can be defined for example as follows: (For $\Omega$, $\partial\Omega=\Gamma_{0}\cup\Gamma_{\boldsymbol{b}}$, $\nu>0$, and $\boldsymbol{f}\in\boldsymbol{L}^{2}(\Omega)$, $\boldsymbol{b}\in\boldsymbol{L}^{2}(\Gamma_{\boldsymbol{b}})$ given) find $\boldsymbol{u}\in\boldsymbol{W}^{1,2}_{\Gamma_{0}}(\Omega)$ and $p\in L^{2}(\Omega)$, such that for all $\boldsymbol{w}\in\boldsymbol{W}^{1,2}_{\Gamma_{0}}(\Omega)$ and $q\in L^{2}(\Omega)$

[TABLE]

Here $\boldsymbol{W}^{1,2}_{\Gamma_{0}}(\Omega)=\{\boldsymbol{v}\in\boldsymbol{W}^{1,2}(\Omega)\,,\;\operatorname{tr}\boldsymbol{v}=\boldsymbol{0}\textrm{ on }\Gamma_{0}\}$ is the standard Sobolev space of functions with $L^{2}$-integrable derivatives in $\Omega$ and the values vanishing on $\Gamma_{0}$. All the data are considered $L^{2}$-integrable here for the sake of brevity. It is worth noting that the above procedure can be used for fluids exhibiting variable viscosity as well; the boundary condition (3) would retain the physical meaning of prescribing the traction on $\Gamma_{\boldsymbol{b}}$.

Since the viscosity is constant, the incompressibility constraint $\operatorname{div}\boldsymbol{u}=0$ allow us to write

[TABLE]

The classical form of the momentum equation (1)1 is thus more usually written as

[TABLE]

Correspondingly, one can write (1) equivalently with $\mathbf{T}$ replaced formally by

[TABLE]

and follow the above procedure, i.e. apply the divergence theorem on $\operatorname{div}\hat{\mathbf{T}}$, the quantity naturally emerging on the boundary being now $-\hat{\mathbf{T}}\boldsymbol{n}$, instead of the traction. Thus, one is led to replace the boundary condition (3) by its analogy

[TABLE]

In the special case that $\hat{\boldsymbol{b}}\equiv\boldsymbol{0}$ this condition is well-known as the do-nothing boundary condition, see [3] for a detailed discussion. In the case that $\hat{\boldsymbol{b}}\neq\boldsymbol{0}$ in general111 By the course of history, the label do-nothing usually refers to the condition (5) with $\hat{\boldsymbol{b}}=\boldsymbol{0}$, but not to the zero traction $\boldsymbol{b}=\boldsymbol{0}$ in (3). The term is not related to the underlying physics in any way. The label originates from the mere fact that the boundary term $(\hat{\boldsymbol{b}},\boldsymbol{w})_{\partial\Omega}$ is absent in the weak formulation (6). Therefore, we do not consider it convenient to use the term do-nothing for the condition (5) with general $\hat{\boldsymbol{b}}$. , we take the liberty to call (5) the full-gradient-traction condition, for the purpose of this paper. The weak solution corresponding to (2) and (5) is then defined as follows: (For $\Omega$, $\partial\Omega=\Gamma_{0}\cup{\Gamma}_{\hat{\boldsymbol{b}}}$, $\nu>0$, and $\boldsymbol{f}\in\boldsymbol{L}^{2}(\Omega)$, $\hat{\boldsymbol{b}}\in\boldsymbol{L}^{2}(\Gamma_{\hat{\boldsymbol{b}}})$ given) find $\boldsymbol{u}\in\boldsymbol{W}^{1,2}_{\Gamma_{0}}(\Omega)$ and $p\in L^{2}(\Omega)$, such that for all $\boldsymbol{w}\in\boldsymbol{W}^{1,2}_{\Gamma_{0}}(\Omega)$ and $q\in L^{2}(\Omega)$

[TABLE]

1.2 Lack of energy estimates, possible remedy and simple flows

The energy estimates for the weak solutions to the steady Navier–Stokes equations are obtained by picking the velocity solution as the test function in the weak momentum equation. By taking $\boldsymbol{w}=\boldsymbol{u}$ in (4) or (6) and using that $\operatorname{div}\boldsymbol{u}=0$ a.e. in $\Omega$, one arrives at

[TABLE]

respectively, allowing for the estimate

[TABLE]

defining either $\tilde{\mathbf{S}}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\mathbf{S}$, $\tilde{\boldsymbol{b}}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\boldsymbol{b}$ and $\Gamma_{\tilde{\boldsymbol{b}}}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\Gamma_{\boldsymbol{b}}$ from (3) and $\sigma\stackrel{{\scriptstyle\mathrm{def}}}{{=}}(2\nu)^{-1}$, or defining $\tilde{\mathbf{S}}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\hat{\mathbf{S}}$, $\tilde{\boldsymbol{b}}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\hat{\boldsymbol{b}}$ and $\Gamma_{\tilde{\boldsymbol{b}}}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\Gamma_{\hat{\boldsymbol{b}}}$ from (5) and $\tilde{\nu}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}{\nu}^{-1}$. Here $\|\cdot\|_{2,\Omega}$ denotes the standard $L^{2}$-norm in $\Omega$ and analogously222We only take $L^{2}$-estimates of the boundary integrals, for simplicity. on $\partial\Omega$. The cubic convective term integral is subject to the following identity, using again that $\operatorname{div}\boldsymbol{u}=0$,

[TABLE]

and is thus determined by the values of velocity on the boundary; its negative represents the inflow rate of the kinetic energy coming through the artificial boundary. In the case of Dirichlet boundary conditions imposed along the entire boundary (or any slip conditions, where $\boldsymbol{u}\cdot\boldsymbol{n}=0$ on $\partial\Omega$), the cubic term would vanish completely and one would obtain an a priori bound for $\tilde{\mathbf{S}}$, concluding the energy estimate for $\boldsymbol{u}$ in $\boldsymbol{W}^{1,2}(\Omega)$ in terms of the data $\boldsymbol{f}$ (by using Poincaré’s or Korn’s inequality, with further requirements on $\Gamma_{0}$). In contrast, the kinetic energy term does not vanish if the artificial boundaries are considered, since the normal component of the velocity is not under control. In this way, for the boundary conditions (3) or (5) with the given boundary data $\tilde{\boldsymbol{b}}\equiv\tilde{\boldsymbol{b}}(\boldsymbol{x})$, one cannot derive the a priori energy bounds.

This fact is reflected by the lack of theoretical results that would guarantee the existence and uniqueness of a steady weak solution. The existence has been stated333Here we only follow the main idea and leave the details on how the results are formulated to the reader. For instance, the referred results are not stated for general $\tilde{\boldsymbol{b}}=\tilde{\boldsymbol{b}}(\boldsymbol{x})$, but rather for the particular choice $\hat{\boldsymbol{b}}=P_{i}\boldsymbol{n}$. only for the case of small data in [3], while without such requirement it remains an open problem in both two and three space dimensions, see also [4, 5]. A detailed discussion can be found in [1], where the uniqueness of a small steady weak solution (for small data) was shown. The uniqueness in the class of all weak solutions (even the uniqueness of the trivial solution for the trivial data) is considered an open problem, as quoted in the first subsection. Similarly, concerning the non-stationary Navier–Stokes problem involving the boundary conditions (3) or (5), the available results only guarantee the existence and uniqueness of weak solutions either for small data or for a small time interval, see [1, 6] and the references therein. In numerical simulations of unsteady flow, this is experienced as a severe instability which appears for larger Reynolds numbers whenever a backflow through such boundary occurs (either through an entire artificial boundary part, or locally, e.g. due to vortices leaving the domain through the outlet) and which may quickly destroy the numerical solution, see also [7].

A number of techniques has been developed to avoid the above mentioned effects in numerical simulations, in particular in the context of unsteady flow, where they are often referred to as “stabilization” methods. We refer the reader to [7] for a selection of various approaches with references. The idea behind one class of these techniques, to illustrate the point briefly, is to modify the boundary conditions (3) or (5) so that the kinetic energy of the eventual inflow (backflow) is compensated. This requires to take the boundary data $\tilde{\boldsymbol{b}}\equiv\tilde{\boldsymbol{b}}(\boldsymbol{x},\boldsymbol{u})$ such that, keeping the notation same as above,

[TABLE]

holds for all admissible functions $\boldsymbol{u}$, say for all $\boldsymbol{u}\in\boldsymbol{L}^{3}(\Gamma_{\tilde{\boldsymbol{b}}})$. In the case of outflow, the cubic term on the right hand side is negative and does not represent any additional restriction on $\tilde{\boldsymbol{b}}$. However, since the occurrence of an inflow ($\boldsymbol{u}\cdot\boldsymbol{n}<0$) cannot be precluded, a suitable quadratic dependence of the boundary data $\tilde{\boldsymbol{b}}$ on the velocity, such as for instance

[TABLE]

is required. In [8], this approach allowed to show the existence of unsteady solutions for arbitrarily large data, and has been utilized in a number of theoretical works since then (see, e.g., Sect. 5 in [9] for a list of further references). This includes the result in [10] showing the existence of a steady weak solution with an arbitrarily large prescribed inflow (imposed via the Dirichlet boundary condition $\boldsymbol{v}=\boldsymbol{g}$, with $\boldsymbol{g}$ given on some part of $\partial\Omega$) subject to the outflow condition

[TABLE]

where $0\leq(\boldsymbol{u}\cdot\boldsymbol{n})^{-}$ denotes the negative part of $\boldsymbol{u}\cdot\boldsymbol{n}$, see also the discussion and the numerical experiments presented in [5]. Analogous results can be shown for fluids with variable viscosity, see for instance [9], where the steady flow of fluids with the pressure- and shear rate-dependent viscosity was considered.

While the above mentioned modified boundary conditions allow to define problems that are mathematically well-posed, it seems that the conditions (3) and (5) are those that are better established and more frequently used in the steady flow numerical simulations, and even for unsteady flows as far as the backflow is not expected or Reynolds number is small. In fact, despite the lack of rigorous results, they both seem quite well-behaved in many problem settings. Note that the full-gradient-traction condition (5) provides a property that all the others mentioned above lack: The condition

[TABLE]

is satisfied by unidirectional simple flows such as Couette or Poiseuille flow, on the inflow and outflow boundaries (defined so that the normal vector $\boldsymbol{n}$ is parallel to the flow). Here $P_{0}\in\mathbb{R}$ is a given constant and, in the resulting simple flow, $p=P_{0}$ on $\Gamma_{\hat{\boldsymbol{b}}}$. As has been discussed in [3] already, the boundary conditions such as (8) result in the velocity solutions with the streamlines distorted near the artificial boundary and the corresponding artefacts in the pressure field. These effects are in many situations undesirable, or are even referred to by some authors as “unphysical”, in particular in the case of flows in channels or pipes, or e.g. in blood vessels, where one naturally supposes that the flow continues undisturbed beyond the artificial boundary. While the “directional” variants such as (9) avoid this disadvantage in the case of outflow, it suffers the same disability in problems, where the inflow boundary needs to be treated in the same way (for one such example, see the plane slider problem studied in [11]).

1.3 The structure of the paper

We present a particular set of flow problems in bounded domains in two or three space dimensions, where one finds multiple steady solutions subject to the same boundary data, the data of the problem being arbitrarily small or even trivial. The specific feature of the presented examples is that the constant traction (3) or full-gradient-traction (5) conditions are prescribed on two444One can experiment with more artificial boundaries as well, see Section 3.6.

boundaries $\Gamma_{i}\subset\partial\Omega$ (one opposite to another), i.e.

[TABLE]

where $\tilde{\boldsymbol{b}}=\boldsymbol{b}$ or $\tilde{\boldsymbol{b}}=\hat{\boldsymbol{b}}$ and where $P_{i}\in\mathbb{R}$ are given constants on each $\Gamma_{i}$. The particularly appealing case where the do-nothing condition $\hat{\boldsymbol{b}}=\boldsymbol{0}$ is imposed on all artificial boundaries is included.

In order to reveal the simple mechanism which promotes this multiplicity, we start in Section 2 by studying the isotropic radial planar flow, where the artificial boundaries are considered at the radii $0<R_{1}<R_{2}$. Given the $P_{i}$–drop555In unidirectional simple flows, such as the flows in channels or pipes, the value $P_{i}$, both whether given as the normal traction or full-gradient-traction, would correspond to the resulting pressure. It is then natural to call the difference $P_{1}-P_{2}$ the pressure drop and to speak about pressure drop problems. However, it is worth noting that in non-unidirectional flows, such as the radial flows, the data $P_{i}$ do not coincide with the resulting pressure. For the sake of brevity, let us denominate $P_{1}-P_{2}$ as $P_{i}$–drop, in this paper. $P_{1}-P_{2}$, the steady problem then reduces to the task of finding one constant, the flow rate $Q\in\mathbb{R}$, which appears to be the solution of a quadratic equation. Apparentally, there is a critical value $P_{\mathrm{crit}}>0$, such that there are two such steady solutions for $P_{1}-P_{2}<P_{\mathrm{crit}}$, which includes the case of trivial boundary data, while for $P_{1}-P_{2}>P_{\mathrm{crit}}$ there is no (isotropic radial) steady solution. Taking the advantage of this reduced setting, we continue in Section 2.3 by studying the unsteady problem (for stationary data). The isotropic radial solution is then given by a single function of time, $Q(t)$, which is found by solving the corresponding ordinary differential equation. We find that the solution either converges asymptotically to one of the steady solutions or it blows up in finite time.

Section 3 is then devoted to the finite element approximation of the flow. We start by examining the isotropic radial setting in Section 3.1, defining the flow domain as one quadrant of a concentric annulus, where the symmetry (or, say, perfect-slip) conditions are imposed on the straight boundaries. We observe that the common simple scheme based on Newton’s method can find both of the two steady solutions, depending on the given initial guess. For large $P_{i}$–drops, a numerical steady solution is found, but it is non-radial. Examining the unsteady case, we confirm numerically the behavior found analytically, including the blow-up of the unsteady solutions in finite time.

Doubt may arise, whether the solutions in the isotropic radial setting are related to the problems where the do-nothing conditions are used most: to the channel-like flows, where the supposed outflow is nearly a unidirectional simple flow. In Section 3.3, we thus make the domain $\Omega$ more narrow and impose no-slip conditions on the walls, studying what resembles a Jeffery–Hamel flow. The numerical results are qualitatively same as in the isotropic radial setting, except that we are unable to find any steady solution for large $P_{i}$–drops.

Ultimately, in Section 3.5 we study the three-dimensional flow through a bifurcating tube, resembling the problems studied in hemodynamics. In this last example, the artificial boundaries are planes placed such that their normal vector coincides with the tube axis (and with the expected flow). In this case, the inflow is imposed by Dirichlet boundary conditions, while the do-nothing condition is prescribed on two “outlet” endings. Analogously to the previous examples, we find two distinct steady solutions, depending on the initial guess assigned to the Newton’s method; one of the solutions displays rapid inflow through the smaller ending.

2 Isotropic radial flow

The multiple solutions described in this paper follow the mechanism which is well illustrated by the following example. The velocity field of the steady isotropic radial planar flow can be written in the form

[TABLE]

The constraint of incompressibility then implies that

[TABLE]

for a constant $Q\in\mathbb{R}$. Here and in what follows we do not indicate that $U$, $U^{\prime}$ are functions of $r$. The convective term takes the form

[TABLE]

Moreover, since

[TABLE]

there holds

[TABLE]

Thus, the balance of momentum in (1) allows us to write the pressure field as

[TABLE]

and is reduced to the following equation for the two unknowns $U\equiv U(r)$ and $P\equiv P(r)$,

[TABLE]

implying that

[TABLE]

Note that the viscous term does not appear here and that the pressure is an increasing function of the radius for any nonzero flux $Q$, no matter of its sign, i.e. for both the inwards or outwards direction of the flow.

Let us define the domain $\Omega$ and its boundary $\partial\Omega=\Gamma_{1}\cup\Gamma_{2}\cup\Gamma_{\mathrm{n}}$ as follows, for given radii $0<R_{1}<R_{2}$ and given $0<\alpha\leq\pi/2$,

[TABLE]

On $\Gamma_{\mathrm{n}}$, the outer normal $\boldsymbol{n}$ is orthogonal to the vector $(x,y)$, implying that $(\nabla\boldsymbol{u})\boldsymbol{n}=(\nabla\boldsymbol{u})^{T}\boldsymbol{n}=U\,\boldsymbol{n}$, so that $-\mathbf{T}\boldsymbol{n}=(P-2\nu U)\,\boldsymbol{n}$. The radial flow (10)–(13) thus meets the perfect-slip boundary conditions (denoting by $\boldsymbol{w}_{\tau}=\boldsymbol{w}-(\boldsymbol{w}\cdot\boldsymbol{n})\boldsymbol{n}$ the tangential part of a vector)

[TABLE]

On $\Gamma_{1}$ and $\Gamma_{2}$, the outer normal $\boldsymbol{n}=(\pm 1/r)\boldsymbol{x}$ is parallel to the direction of the flow. Therefore, $(\nabla\boldsymbol{u})\boldsymbol{n}=(\nabla\boldsymbol{u})^{T}\boldsymbol{n}=(U+rU^{\prime})\boldsymbol{n}$ and the observed traction is given by

[TABLE]

In particular, for the incompressible fluid, (10) and (14) lead to the following relation for the normal traction,

[TABLE]

where $r=R_{i}$ on $\Gamma_{i}$, $i=1,2$.

2.1 Nontrivial steady solution for trivial data

As a consequence of (18), the incompressible Navier–Stokes equations in $\Omega$ subject to the perfect-slip boundary conditions (16) and the zero traction boundary conditions (3)

[TABLE]

exhibit two isotropic radial solutions: one trivial, defined by $P_{\infty}=Q=0$, and one given by

[TABLE]

In the case that the zero traction conditions are replaced by the do-nothing conditions (5)

[TABLE]

since by simple computation $p\boldsymbol{n}-\nu\,(\nabla\boldsymbol{u})\boldsymbol{n}=-\mathbf{T}\boldsymbol{n}+\nu\,(\nabla\boldsymbol{u})^{T}\boldsymbol{n}=(P_{\infty}-\tfrac{1}{2}r^{-2}Q^{2}+\nu r^{-2}Q)\,\boldsymbol{n}$, one makes the similar observation, the non-trivial solution being defined by

[TABLE]

The resulting flux observed for the do-nothing conditions is half of the flux for the zero traction conditions, cf. (7).

2.2 Two, one or no solution for given $P_{i}$–drop

Let us further consider the $P_{i}$–drop problem666By the phrase $P_{i}$–drop we avoid using the more common, but in our setting confusing, phrase pressure drop. See the footnote in Section 1.3. , where the flow is subject to two values of the normal traction prescribed on $\Gamma_{1}$ and $\Gamma_{2}$,

[TABLE]

For the given radii $R_{i}$ and the normal tractions $P_{i}$, one obtains by subtracting (18) for $i=1,2$, the following quadratic equation for the unknown $Q$,

[TABLE]

Having $Q$ in hand, $P_{\infty}$ is computed by using (18) for either $i=1$ or $2$.

For $P_{1}-P_{2}<P_{\mathrm{crit}}$, there are two solutions for $Q$ as follows,

[TABLE]

It is worth noting that the larger of the two solutions $Q$ is always positive and that it is a decreasing function of the difference $P_{1}-P_{2}$. For the special case $P_{1}-P_{2}=P_{\mathrm{crit}}$, there is only one solution, given by $Q=2\nu$. Interestingly, there is no steady isotropic radial solution once the $P_{i}$–drop is too large, $P_{1}-P_{2}>P_{\mathrm{crit}}$.

In the case of full-gradient-traction conditions

[TABLE]

one observes the analogous result, the critical $P_{i}$–drop being then at

[TABLE]

and the unique solution being then given by $Q=\nu$, for $P_{1}-P_{2}=\hat{P}_{\mathrm{crit}}$.

2.3 Stability of the one-dimensional problem

Let us now consider the unsteady flow of the same type, i.e. let us find

[TABLE]

that solve the Navier–Stokes system of equations

[TABLE]

supplemented by the (constant in time) boundary conditions (16) and (19) or (21). Same as in the stationary case, the incompressibility constraint implies

[TABLE]

and, by using (11) and (12) the momentum equation implies that

[TABLE]

which gives the pressure field

[TABLE]

Applying the constant $P_{i}$–drop by imposing (19), one obtains from (17) the following two equations for the two unknown functions $P_{\infty}(t)$ and $Q(t)$,

[TABLE]

By subtracting, the following single equation for $Q(t)$ is concluded,

[TABLE]

with $P_{\mathrm{crit}}$ from (20). Its solution $Q(t)$ appears in three different forms, depending on the sign of $c$.

For $c>0$, i.e. for $P_{1}-P_{2}<P_{\mathrm{crit}}$, the two stationary solutions $Q=b\pm\sqrt{c}$ are indeed found same as in the previous subsection, while all other solutions are of the form

[TABLE]

with $t_{0}\in\mathbb{R}$ to be determined by the initial condition. In particular, we conclude the following asymptotics of the maximal solutions starting from $Q(0)$,

[TABLE]

We have found that one of the stationary solutions is asymptotically stable within the considered class of radial solutions, while the second one is unstable. In particular, in the case of zero $P_{i}$–drop, the obvious trivial solution $Q\equiv 0$ is asymptotically stable in the above sense, while the other solution given by $Q\equiv 4\nu$ is unstable.

In the special case $c=0$, i.e. for $P_{1}-P_{2}=P_{\mathrm{crit}}$, there is the single stationary solution $Q\equiv 2\nu$, while the other solutions are

[TABLE]

with $t_{0}\in\mathbb{R}$. In this case,

[TABLE]

again with $0<t_{0}<+\infty$ depending on $Q(0)$.

Finally, for $c<0$ ($P_{1}-P_{2}>P_{\mathrm{crit}}$), where no stationary isotropic radial solution can be found, all solutions are in the form

[TABLE]

For any initial value, $Q(t)$ will blow up to $+\infty$ in finite time (obviously less than $\tfrac{\pi a}{\sqrt{-c}}$).

Same as before, in the case of the full-gradient-traction $P_{i}$–drop (21), the analogous observations hold with the only difference that $b$ and $c$ in (24) are to be replaced by

[TABLE]

with $\hat{P}_{\mathrm{crit}}$ from (22). We note in particular the case of do-nothing conditions on both boundaries, i.e. $P_{1}=P_{2}=0$, where the flow decays to zero if $Q(0)<2\nu$, while it blows up in finite time for any initial $Q(0)>2\nu$.

3 Numerical simulations

We performed a number of numerical simulations of the flows subject to the constant traction or full-gradient-traction conditions on opposite artificial boundaries. In what follows we present some of them, on some we report only briefly, so that the presentation is confined within a reasonable extent. The numerical results were obtained by the finite element scheme implemented using FEniCS [12, 13]. Navier–Stokes equations were discretized using Taylor–Hood finite element pair on a triangular mesh within the mixed formulation, i.e. we solve for $(\boldsymbol{u}_{h},p_{h})\in V_{h}\times Q_{h}$ such that (taking $\boldsymbol{f}=\boldsymbol{0}$, for simplicity)

[TABLE]

holds for all $(\boldsymbol{w},q)\in V_{h}\times Q_{h}$, cf. (4) and (6). Either the constant traction (19) is prescribed, $\tilde{\boldsymbol{b}}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\boldsymbol{b}=P_{i}\boldsymbol{n}$ on $\Gamma_{i}$, in which case $\tilde{\mathbf{S}}(\nabla\boldsymbol{u}_{h})\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\nu\,(\nabla\boldsymbol{u}_{h}+(\nabla\boldsymbol{u}_{h})^{T})$, or the weak formulation with the full gradient is used, $\tilde{\mathbf{S}}(\nabla\boldsymbol{u}_{h})\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\nu\nabla\boldsymbol{u}_{h}$, and the corresponding conditions (21) are imposed, prescribing again $\tilde{\boldsymbol{b}}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\hat{\boldsymbol{b}}=P_{i}\boldsymbol{n}$ on $\Gamma_{i}$.

The Taylor–Hood finite element function spaces $V_{h}\subset W^{1,2}(\Omega)^{2}$ and $Q_{h}\subset L^{2}(\Omega)$ are defined in a standard way, see e.g. [13]. On the remaining parts of the boundary we prescribe either (in Sections 3.1 and 3.2) the perfect-slip boundary conditions (16), in which case $V_{h}\subset\{\boldsymbol{v}\in W^{1,2}(\Omega)^{2}\,,\;\operatorname{tr}(\boldsymbol{v}\cdot\boldsymbol{n})=0\textrm{ on }\Gamma_{\mathrm{n}}\}$; or (in Sections 3.3–3.5) the no-slip condition on $\Gamma_{0}\subset\partial\Omega$, taking $V_{h}\subset\{\boldsymbol{v}\in W^{1,2}(\Omega)^{2}\,,\;\operatorname{tr}\,\boldsymbol{v}=\boldsymbol{0}\textrm{ on }\Gamma_{0}\}$.

While focusing on steady problems, we present a few unsteady simulations as well. The unsteady problems (in Sections 3.2 and 3.4) are first discretized by the fully implicit (backward Euler) scheme with the fixed time step $\tau>0$, i.e. we solve subsequently at each time level $t_{k}=k\tau$, $k=1,2,\ldots$, the following equation for the unknowns $(\boldsymbol{u}_{h}^{k},p_{h}^{k})\equiv(\boldsymbol{u}_{h}(t_{k}),p_{h}(t_{k}))\in V_{h}\times Q_{h}$,

[TABLE]

for all $(\boldsymbol{w},q)\in V_{h}\times Q_{h}$, where $\boldsymbol{u}_{h}^{0}=\boldsymbol{u}_{h}(0)$ is a given initial condition. Both the steady and the unsteady problem result in (one, or a sequence of) nonlinear algebraic systems, that are solved by the FEniCS built-in DOLFIN Newton solver [13, 14] with the exact Jacobian matrix computed by automated differentiation. All numerical results presented were checked to be independent of a further refinement in $h$ or $\tau$.

3.1 Steady isotropic radial setting

We start by performing numerical simulations in the setting corresponding to the steady isotropic radial flow discussed in Section 2.2. We define the flow domain $\Omega$ as in (15) with $\alpha=\pi/2$, so that the constraint $\boldsymbol{u}\cdot\boldsymbol{n}=0$ on $\Gamma_{\mathrm{n}}$ can be written, denoting $\boldsymbol{u}=(u,v)$, as

[TABLE]

which simplifies the implementation. We set $R_{1}=1$, $R_{2}=3$ and $\nu=1$.

Prescribing the constant traction (19), $\boldsymbol{b}=P_{i}\boldsymbol{n}$ on $\Gamma_{i}$, for a set of $P_{i}$–drops $P_{1}-P_{2}$, we report the resulting normal flux777 Note that $\boldsymbol{e}_{r}=\frac{1}{r}(x,y)=-\boldsymbol{n}$ on $\Gamma_{1}$ while $\boldsymbol{e}_{r}=\boldsymbol{n}$ on $\Gamma_{2}$, and that the definition of $Q_{\mathrm{num}}$ is independent on whether $i=1$ or $2$, due to $\operatorname{div}\boldsymbol{u}=0$. For the isotropic radial flow presented in Section 2, where $\boldsymbol{u}=Q\tfrac{1}{r}\boldsymbol{e}_{r}$, there is $Q_{\mathrm{num}}=Q$.

[TABLE]

in Figure 1a. We make the following observation.

• •

For $-1<P_{1}-P_{2}\leq 1.\overline{7}$, by starting Newton’s method from the zero initial guess we obtain the isotropic radial solutions, recovering the lower branch of (20). Note that the critical $P_{i}$–drop with only one corresponding radial solution derived in Section 2.2 is $P_{\mathrm{crit}}=2\nu^{2}\,(R_{1}^{-2}-R_{2}^{-2})=16/9=1.\overline{7}$.

• •

Starting Newton’s method from the velocity field with large enough radial flux, specifically starting it from the isotropic radial solution $\boldsymbol{u}=r^{-2}Q_{\mathrm{init}}\boldsymbol{x}$ with $Q_{\mathrm{init}}=3$, for the same range of $P_{i}$–drops as above we obtain different numerical solutions, recovering the upper branch of (20).

• •

For $P_{1}-P_{2}$ larger than but close to the critical value $1.\overline{7}$, Newton’s iterations (started by, for instance, $Q_{\mathrm{init}}=0$ or $Q_{\mathrm{init}}=3$) get lost. The algebraic residuum oscillates and does not show much hope for convergence.

• •

For the $P_{i}$–drop given large enough, say $P_{1}-P_{2}\geq 2$, in which case there exists no isotropic radial solution, we are finding non-radial numerical solutions. For convenience, we present a simple “non-radiality indicator”

[TABLE]

in Figure 1b. By naive numerical continuation, namely by decreasing the $P_{i}$–drop subsequently while using the previous non-radial solution as the initial guess, we are finding the non-radial solutions for $P_{1}-P_{2}$ down to $1.65$, see Figure 1a. Note on Figure 1b that the non-radiality gradually vanishes with decreasing $P_{i}$–drop, the solutions eventually joining the lower branch of isotropic radial solutions.

Numerical simulations using the full-gradient-traction (21), $\hat{\boldsymbol{b}}=P_{i}\boldsymbol{n}$ on $\Gamma_{i}$, lead to analogous results, as reported in Figure 1c and 1d.

• •

Note that the critical $P_{i}$–drop here equals $P_{\mathrm{crit}}=\tfrac{1}{2}\nu^{2}\,(R_{1}^{-2}-R_{2}^{-1})=4/9=0.\overline{4}$.

• •

The main difference regards the non-radial solutions found (e.g., first from the zero initial guess for $P_{1}-P_{2}=2$, and then numerically continued to lower $P_{i}$–drops). As shown in Figure 1c, the non-radial solutions now live down to the negative $P_{1}-P_{2}<-1.3$.

• •

In particular, we thus found three888Or four, actually: the non-radial solutions appear in two variants symmetric around the axis $x=y$. distinct numerical solutions for trivial data, i.e. for the do-nothing boundary condition prescribed on $\Gamma_{0}$ and $\Gamma_{1}$. The velocity and pressure fields of the corresponding non-radial solution are shown in Figure 2.

3.2 Unsteady isotropic radial setting

Similarly as in Section 2.3, we are interested in the unsteady flows in the very same setting as above, i.e. we keep the geometry, the boundary conditions and the viscosity same as in the previous section, but we look for the unsteady solutions to (23) with the given initial state

[TABLE]

We only perform simulations with stationary boundary conditions.

First, we set the initial velocity field simply as the isotropic radial velocity

[TABLE]

Running the numerical simulations for the $P_{i}$–drops $P_{1}-P_{2}=-1$, [math],…,$3$ and $Q_{0}=-1$, [math],…,$5$, we always reproduced the exact unsteady isotropic radial solutions derived in Section 2.3. Note that we stop the simulations at $t=10$, or if the velocity is becoming large, namely as soon as $Q_{\mathrm{num}}\geq 20$.

In particular, the numerical solutions remain isotropic radial and show the blow-up in finite time for $P_{1}-P_{2}>P_{\mathrm{crit}}$; while for $P_{1}-P_{2}<P_{\mathrm{crit}}$, where there exist two steady solutions, the numerical solutions either converge asymptotically to the smaller (stable) steady isotropic radial solution, or they blow up in finite time, depending on the initial state.

3.3 Steady $P_{i}$–drop diverging channel flow

Inspired by the behavior of isotropic radial flows, we now focus on the simulations which are not restricted to such simplified setting. We start, for simplicity, by keeping the domain geometry $\Omega$ similar as before, only with $\alpha=\pi/8$ in (15), while we keep $R_{1}=1$ and $R_{2}=3$. We now prescribe the no-slip boundary conditions on the straight walls,

[TABLE]

We set the viscosity $\nu=1/20$, so that the resulting fluxes of the two branches of solutions will be of the same order as in Section 3.1. We report them in Figure 3a, showing the behavior analogous to the one presented in Figure 1. The results for traction and full-gradient-traction boundary data are now plotted together, showing that

• •

again, the initial guess based on $Q_{\mathrm{init}}=3$ reveals the upper branch of steady solutions;

• •

the resulting flux on the upper branch for given traction is larger than the one for full-gradient-traction, although the difference is less pronounced than it was for the isotropic radial flow;

• •

no numerical steady solutions are obtained for large enough $P_{i}$–drops.

For illustration, we present in Figure 4 the velocity and pressure fields of the non-trivial steady numerical solution for zero full-gradient-traction $P_{i}$–drop, i.e. with the do-nothing boundary conditions. It shows that the pressure is nearly zero in the vicinity of the artificial boundaries, as expected, while it is lower inside the domain. The pressure field indicates that the velocity cannot be radial. The flux is around $Q_{\mathrm{num}}=2.7$, in accordance with Figure 3a.

3.4 Unsteady $P_{i}$–drop diverging channel solutions

Analogously to Section 3.2, we look for the unsteady solutions, while keeping the geometry, boundary conditions and viscosity the same as in Section 3.3. For simplicity, we set the initial velocity field numerically as the isotropic radial velocity (27) (despite the fact that it does not satisfy the no-slip boundary conditions). Running the numerical simulations for $Q_{0}=-1$, [math], …, $4$ and various $P_{i}$–drops, we obtain results analogous to those in Section 3.2. We report in Figure 3b the solutions for three full-gradient-traction $P_{i}$–drops. For $P_{1}-P_{2}=0$, all but one solutions converge to zero, while the one starting from $Q_{0}=4$ displays the blow-up in finite time. Similarly for $P_{1}-P_{2}=20$, the solutions either converge to the lower steady solution, or blow up if the initial-state flux is larger than that of the upper steady solution. For $P_{1}-P_{2}=30$, all the solutions blow up, the fluxes $Q_{\mathrm{num}}(t)$ appear to be only shifted in time.

3.5 Flow through a bifurcating tube

Finally, we study the three-dimensional flow in the geometry of a bifurcating tube, let us say a bifurcating blood vessel, see Figure 5. The inflow through one inlet boundary (on the top) is given, by prescribing the parabolic velocity profile via Dirichlet boundary condition, the no-slip conditions are assumed on the walls. The outflow is possible through two outlets (on the bottom), where we set the do-nothing boundary condition. The geometry of the domain was constructed such that the artificial boundaries represent circular cross-sections of the vessel parts, the diameters of the inlet and the two outlets being $0.5$, $0.2$ and $1.0$. The length of the domain (from the inlet to the outlets) is $3$ and the distance of the outlets centers is approximately $1$, the inflow velocity is $1$ at its midflow maximum and the viscosity999The quantities we work with can be considered as non-dimensional here, the viscosity thus represents a reciprocal of Reynolds number. Note that Reynolds number $\mathrm{Re}=50$ is related to the characteristic length equal to the diameter of the larger outlet, or to the distance of the centers of the two outlets, and to the characteristic velocity equal to the peak velocity at the inlet. Note that the peak backflow velocity of the second solution reaches five times more, which would then correspond to $\mathrm{Re}=250$. is given by $1/\nu=50$.

The fact that the cross-sectional area of the two outlets differ significantly is important for the existence of two steady laminar numerical solutions. Figures 5a and 5b present two steady flows subject to the same boundary conditions: the no-slip condition on the walls and the do-nothing on the outlets. The non-trivial flux through the inlet was prescribed only in order to make the example more representative for the real application; prescribing the trivial data on the inlet as well would result in the very same non-unique behaviour.

3.6 Comment on other examples

Analogously, one can obtain more than two solutions in domains with more than two artificial boundaries treated with the conditions (3) or (5), provided that the ratios of their area are large enough and the domain allows for flows with small enough dissipation. One such example, which we do not report in detail here, would consist of a triangular domain, where (3) or (5) is given on certain parts (of a different length) of each edge of the triangle. For certain lengths of the artificial boundaries, it is easy to find combinations of $P_{i}$–data allowing for up to three (reasonably small, laminar) steady solutions. Such a setting, however, is not of much practical importance.

On the other hand, by for example gradually prolonging the outlets in Figure 5 by straight tubes, while watching the resulting inflow through the smaller outlet of the solution illustrated in Figure 5b, one would observe a significant increase of the inflow velocity, loosing soon the possibility of finding two (laminar) steady solutions. This is in accordance with the fact, that the do-nothing conditions are used in practical simulations on regular basis without reporting multiple solutions such as those considered in this paper. According to our experience, it is only in specific applications that one can encounter this non-uniqueness accidentally.

4 Conclusion

The fact that the boundary conditions prescribing a constant value of traction or what we call the full-gradient-traction here, including the do-nothing boundary conditions, inhibit from establishing the well-posedness theory for the steady flow problem, is well known by the mathematical community concerned with the flows of incompressible fluids. The banality with which the do-nothing conditions are used in a multitude of numerical simulations none the less, reflects the amenity of these boundary conditions and their robustness in many practical flow problems. It is not the aim of this paper whatsoever to reason against using them, rather we attempt to contribute to the discussion on how to use them and what results can be expected.

Clearly, following the examples presented in this paper, the uniqueness of the steady solution to the problems subject to boundary conditions that include (3) or (5) cannot be expected in general, not even for small or trivial data. Although this may not be surprising, it has not been reported in the literature in this detail before, to the best of our knowledge. The more prominent question, whether multiple solutions can be found in problems where only one connected part of the boundary is treated in this way, remains unresolved. We can only report without any details that we have tried to find such numerical examples as well, yet with no success.

Inspired by Sections 3.2 and 3.4, other questions arise as well: in the case of multiple steady solutions (say, for small data), is there still a unique stable steady solution, as both the numerical simulations in the diverging channel and the analytical isotropic radial solutions suggest? In particular, is the unique small solution established by Galdi in [1] the one (and only one) that is stable? Somewhat related is the question, which of the multiple solutions can be delivered by various non-linear algebraic problem solvers, as we only reported the numerical experiments with Newton’s method (which can reveal multiple solutions, depending on the initial guess), and the solutions to the unsteady flow (which, if considered as a tool for finding the steady solution, either failed or revealed one steady solution at most).

Other presented examples suggest the lack of existence of the steady solution for certain problems with large enough boundary data. Such behavior of Navier–Stokes equations is expected in general, but largely due to the onset of an unstable non-laminar flow regime, in which case the existence of the numerical solution can be affected by various convective term stabilization techniques. Here, the isotropic radial example reveals a rather different mechanism related to the blow-up of the unsteady flow in finite time due to the kinetic energy being brought into the domain by the inflow through one of the artificial boundaries. It seems obvious that this behavior can only be avoided by altering the boundary conditions.