Variational resolution of outflow boundary conditions for incompressible Navier-Stokes

This paper focuses on the so-called Weighted Inertia-Dissipation-Energy (WIDE) variational approach for the approximation of unsteady Leray-Hopf solutions of the incompressible Navier-Stokes system. Initiated in [56], this variational method is here extended to the case of non-Newtonian fluids with power-law index $r\geq11/5$ in three space dimensions and large nonhomogeneous data. Moreover, boundary conditions are not imposed on some parts of boundaries, representing, e.g., outflows. Correspondingly, natural boundary conditions arise from the minimization. In particular, at walls we recover boundary conditions of Navier-slip type. At outflows and inflows, we obtain the condition $-\frac12|\boldsymbol{v}|^2\boldsymbol{n}+\mathbb{T}\boldsymbol{n}=0$. This provides the first theoretical explanation for the onset of such boundary conditions.


Introduction
In this work, we are interested in a variational resolution technique of the incompressible Navier-Stokes system by means of the so-called Weighted Inertia-Dissipation-Energy (WIDE) functional approach [63,66].
We consider the unsteady flow of an incompressible fluid through a generalized channel Ω with inlets Γ D , walls Γ N , and outlets Γ i F , i = 0, . . . , n, as depicted on Figure 1. The flow is modeled in the bulk by the incompressible Navier-Stokes system div v = 0, ∂ t v + v · ∇v − div S(Dv) + ∇p = f , (1.1) where S describes the constitutive relation between the Cauchy stress tensor and the symmetric velocity gradient Dv := 1 2 (∇v + (∇v) T ). The main focus of this paper is to show that weak solutions to this system can be obtained as limits as ε → 0 + of minimizers of the WIDE functionals (1.2) when minimized over whole trajectories satisfying only the constraints div v = 0, v| ΓD = v D , v| ΓN · n = 0, for some given (inflow) data v D and some given net flux rates F i through Γ i F , i = 0, 1, . . . , n. The symbol • v denotes a kind of material derivative defined by where ε ijk is the Levi-Civita symbol. Further, the symbol n is the outward normal vector on ∂Ω. In fact, • u can be replaced by the usual material derivative • u := ∂ t u+u·∇u in some specific cases, as explained below. Moreover, the function s is an analogue of S on the boundary, modeling the friction forces on Γ N . The functions S ε and s ε approximate S and s, respectively, by improving their asymptotic growth.
The relation between the minimization of the WIDE functional I ε and the Navier-Stokes system is revealed by formally computing the Euler-Lagrange equation of I ε , that is, by assuming smoothness. Let v ε minimize I ε in some reasonable set of trajectories obeying (1.3) and compute the Gateaux derivative of I ε with respect to v in a direction ϕ, which does not have fixed boundary values on Γ N and Γ F . Via integration by parts (more precisely, the Stokes theorem) this leads to (1.5) where q ε is the Lagrange multiplier corresponding to the constraint div v ε = 0. Here, the subscript τ denotes the tangential part of a vector on ∂Ω, i.e., w τ := w − (w · n)n = n × w × n.
As seen, the minimization of I ε corresponds to an elliptic-in-time regularization of the incompressible Navier-Stokes system. Such regularizations in the setting of linear and nonlinear parabolic systems are quite classical and have already been considered in [44,33,55], especially as tool for tackling regularity issues. The reader is referred to the classical monograph [47] for an account of results in the linear setting. The Navier-Stokes system has also been investigated by elliptic-intime regularization [45,46], in a setting however which does not admit a variational structure.
The application of the WIDE variational approach to incompressible Navier-Stokes has been initiated in [56]. There, no-slip boundary conditions are imposed on the whole boundary ∂Ω and the fluid is assumed to be Newtonian. The present paper extends the reach of [56] by allowing for in-and outlets, by considering general, nonhomogeneous boundary conditions, and by allowing non-Newtonian effects in the fluid.
The possibility of dealing with nonhomogeneous boundary conditions, especially with the outflow boundary condition (1.11), is very relevant in relation with applications, see [29] and references therein. To the best of our knowledge, conditions (1.11) are however still missing a thorough physical justification. An important contribution of this paper is hence that of providing a variational justification of boundary conditions (1.11), for they arise as natural conditions via minimization of I ε .
The scope of the paper is to make the above formal argument rigorous: under suitable assumptions on boundary and initial data, we prove that the minimizers of I ε over trajectories constrained by (1.3) converge weakly for ε → 0 + (up to subsequences) to a Leray-Hopf solution [67] of the Navier-Stokes system with boundary conditions (1.9)-(1.11). As mentioned, the accent is here not on existence, which for this system is already known, see [15,34], but on variational charaterization of natural boundary conditions. The structure of the paper is as follows. We collect some detail on the physical setting of the problem and on the WIDE functional approach in Section 2. Notation and assumptions are then presented in Section 3. The statement of our main result, Theorem 1, is in Section 4, where we also record a collection of remarks in order to put the statement in context. Section 5 eventually contains the proof of Theorem 1.

Physical motivation
Before we proceed with the rigorous mathematical treatment, let us elaborate on the physical meaning of the whole procedure and its relation with the existing theory.
In Figure 1, the set Ω is a representative of an open bounded subset of R d , d ≥ 2, whose Lipschitz boundary ∂Ω is divided into three different types: Γ D : A Dirichlet boundary condition is prescribed here. This can model an adhesive boundary (no-slip), or (more importantly) it represents a prescribed inflow in the nonhomogeneous case. Γ N : This set represents the impermeable walls, where v · n = 0. No condition is imposed in the tangential direction a priori. Γ i F : The fluid passes freely through this boundary, only the corresponding net flux must be equal to a given amount, expressed by a function F i of time. As such, this represents an artificial boundary, such as outlet or inlet, where nothing is known a priori about the flow, except for the net flux. Note that due to the incompressibility, it is enough to prescribe F i for i = 1, . . . , n only, for example.
The sets Γ D , Γ N , and Γ i F , i = 0, 1, . . . , n, are open, but possibly not connected (in order to allow multiple walls, inflows/outflows with shared flux rate etc.). However, for simplicity we shall assume that they consist of a finite number of connected components. Of course, it is possible to simplify this very general setting by omitting certain types of boundaries. We will nevertheless assume that the boundary setting allows for a corresponding Poincaré inequality in Ω.
2.1. The problem of outflow boundary conditions. Finding a reasonable outflow boundary condition for a flow of an incompressible fluid through a channel is a longstanding problem, which is obviously of great importance in numerics and applications. On one hand, some outflow boundary condition seems to be required in order to make the problem well-posed. On the other hand, usually no a priori information is available at the outflow. Consider for example, the flow in a section of a pipe. One would formulate the outflow boundary condition by trying to replicate the flow behavior as if the pipe was a section of a much longer pipe, an idealization which is often interpreted as the size of the system scaled up. To this day, there is still no satisfactory theory dealing with this problem. Most of the existing mathematical theories of (incompressible) fluids consider internal flows only, an assumption that is actually very rarely met in reality. This is in a sharp contrast with the numerical simulations of the incompressible flow, where several types of outflow boundary conditions are used.
Starting from [26], where the so-called do-nothing boundary condition − pn + ν(∇v)n = 0 on Γ 0 F , ν > 0, (2.1) was introduced for the first time, many outflow boundary conditions have been proposed. All these can be roughly divided into two classes: (A) Outflow boundary conditions that are made to fit the geometry of the problem (experiment) at hand. This leads to no versatility of application and usually also to poor mathematical properties. (B) Outflow boundary conditions that are guessed, based on the required mathematical properties (such as the validity of the energy estimate, controllable backward flow etc.).
Although the do-nothing boundary condition (2.1) is often deemed to be natural (since it eliminates the whole boundary term in the variational formulation of the Navier-Stokes equations), this condition is actually a canonical representative of the class (A). Indeed, this condition is compatible with a Poiseuille flow through a straight channel. This is however the effect of the term (∇v)n, which vanishes if the outlet is perpendicular to the direction of the flow, which is clearly a very geometry-dependent property. See [29] for some prototypical examples, where (2.1) leads to unphysical flows. We also remark that no energy estimate for standard weak solutions of the Navier-Stokes system with (2.1) imposed on a part of the boundary is available, and hence no corresponding large-data, global-in-time existence theory.
On the other hand, the outflow boundary conditions of the class (B) allow to control the unsigned term |v| 2 v · n on the outlet, hence giving an energy estimate. The outflow boundary conditions studied, e.g., in [29], [12], or [13], are typical representatives of the class (B). These are also particular cases of the class of energy-preserving boundary conditions discovered and studied in [15] or more recently in [54] using a different approach. In these works it is also shown that some of those boundary conditions may produce reliable numerical results with respect to experiments, even in the case of a turbulent flow. Unfortunately, the approach based on [15] or [54] provides no physical motivation for such boundary conditions. This indeed presents a problem since the class (B) is actually huge (basically due to the incompressibility constraint, cf. [35]) and different outflow boundary conditions can lead to dramatically different flow characteristics, at least near the outlet. What is seemingly missing is a specific selection criterion. 2.
2. An optimization problem in order to qualify boundary conditions. In this work, we attack the problem of outflow boundary conditions by a rather unorthodox method, that has been so far used only in [9] to derive the boundary conditions of the do-nothing type for the stationary Stokes system. Here, we model the flow by the more appropriate unsteady Navier-Stokes equations and we also consider that the fluid is non-Newtonian (for example, we allow the dependence of the viscosity on the shear rate). The underlying mathematical idea is that variational (weak) formulations of PDEs (which are also physically more natural) may implicitly encode boundary conditions if the test functions have enough freedom near the boundary. This leads to looking for a solution in some large space with unspecified boundary conditions. It thus seems natural to reformulate everything as an optimization problem in such a space. Choosing an appropriate functional I to minimize is, of course, a non-trivial task. On the other hand, in case of the dissipative systems such as the Navier-Stokes equations, the energy dissipation itself may serve as a good candidate. In fact, this provides a possible physical explanation and a selection criterion for the obtained (outflow) boundary conditions: they are such that the corresponding flow dissipates the least amount of energy, in some sense. This approach has several other advantages. Firstly, if the functional I is indeed related to the physical energy, then one readily gets an energy estimate and, consequently, the existence of a (weak) solution follows. This is not in contradiction with the lack of an energy estimate for (2.1) since we are able to derive (2.1) by our method only for certain shear-thickening fluids, cf. (4.18) below. Another motivation for prescribing the outflow boundary conditions in such an implicit way can be found in [57]. There, it is argued that just by applying a finite-element discretization to the variational formulation of the system with the so-called no boundary condition, one implicitly prescribes some outflow boundary conditions that turn out to possess superior numerical properties. See also [27] and [60] for an explicit resolution of such boundary conditions in 1D.
Adapting the idea of [9] to our setting is not a trivial task. Indeed, we need to tackle the additional problems of adding the time evolution and of the intrinsic nonlinearity of the Navier-Stokes equation (which cannot be easily treated as a constraint). It turns out that both these issues can be solved by a clever choice of the functional I, whose Euler-Lagrange equations coincide (in the bulk) with the unsteady Navier-Stokes system in the sense of a certain limit.
2.3. The WIDE functional. As already remarked in the Introduction, the Euler-Lagrange equation (1.5) of the WIDE functional I ε is nothing but an elliptic-intime regularization of the original Navier-Stokes system (1.1). The use of WIDE variational approach can be tracked back at least to Ilmanen [31], who used it to tackle existence and partial regularity of the Brakke mean-curvature flow of varifolds. An application to existence of periodic solutions for gradient flows is given by Hirano [30]. The variational nature of elliptic regularization is at the core of [23, Problem 3, p. 487] of the classical textbook by Evans. Two examples of relaxation related to micro-structure evolution have been provided in [17] and the case of mean-curvature evolution of Cartesian surfaces is in [65]. The analysis of the WIDE approach for abstract gradient flows for λ-convex and nonconvex energies is in [52,6] in the Hilbertian case and in [61,62] in the metric case. Melchionna [48] extended the theory to classes of nonpotential perturbations and Bögelein, Duzaar, & Marcellini [10] used this variational approach to prove the existence of variational solutions to the equation u t − ∇ · f (x, u, ∇u) + ∂ u f (x, u, ∇u) = 0 where the field f is convex in (u, ∇u).
Doubly nonlinear parabolic evolution equations have been tackled by the WIDE variational formalism as well. The first result in this direction is by Mielke & Ortiz [50], where the case of rate-independent processes is addressed. The corresponding time discretization has been presented in the subsequent [51] and an application to crack-front propagation in brittle materials is in [39]. The rate-dependent case has been analyzed in [2,3,4,5]. See also [41] for a stability result via Γ-convergence [18] and [49] for an application to the study of symmetries of solutions.
In the dynamic case, De Giorgi conjectured in [21] that the WIDE functional procedure could be implemented in the setting of semilinear waves. This has been ascertained in [66] (for the finite-time case) and by Serra & Tilli [63] (for the infinite-time case). The possibility of following this same variational approach in other hyperbolic situations has also been pointed out in [21]. Indeed, extensions to mixed hyperbolic-parabolic semilinear equations [42], to different classes of nonlinear energies [43,64], and to nonhomogeneous equations [68,69] are also available. The validity of the WIDE approach to the wave equation in time-dependent domain [19], dynamic perfect plasticity [20]. Eventually, the incompressible Navier-Stokes system has been tackled in [56].
We build on the work [56], where it is shown that the Navier-Stokes equations can be obtained as limit of Euler-Lagrange equations of mimimizers of certain WIDE functionals, constrained to v = 0 on the boundary. Here, we remove this constraint and, in effect, we leave to the functional to decide which boundary condition is optimal. Moreover, we also allow the non-Newtonian effects both in the bulk and on the boundary Γ N . This leads to the functional I ε defined in (1.2). The positive number ε tends to zero and it parametrizes the weight e − t ε , which dampens the impact of the future evolution on a current state of the given system. The function f represents an external body force density, such as gravity. Further, the function S : R 3×3 sym → R 3×3 sym represents the stress-strain relation and the function S ε is given by S ε (A) := S(A) + εσ 4 |A| 2 A + εσ q |A| q−2 A, σ 4 , σ q > 0, q > 1. (2. 2) The corresponding integral in the definition of I ε is the amount of the dissipated energy due to internal friction. We postpone the formulation of precise assumptions on S to the next section, but for the time being, the reader may think of, e.g., the Ladyzhenskaya model satisfying obviously which is actually a canonical representative of the class of models that is considered later. In the special case r = 2, it is assumed that S is linear, i.e., S(A) = 2σ 2 A corresponding to the classical Navier-Stokes model. The nonlinear model (2.3) is plausible also if 1 < r < 2, leading to shear thinning fluids, but this case is excluded in the analysis below. The presence of ε > 0 in S ε further improves its asymptotic growth. This helps us in several places to deal with the convective term in the equation. The parameters σ 2 , σ r , σ 4 , and σ q can be interpreted as of the functions s and s ε is completely analogous to that of S and S ε , respectively, and we shall impose s ε (u) := s(u) + ερ 4 |u| 2 u + ερ q |u| q−2 u, ρ 4 , ρ q > 0, q > 1, (2.4) where, for the present time, we choose We do not claim that I ε is easy to justify physically. At least some of its terms however bear a clear physical meaning. First of all, the term −f · v obviously corresponds to the work done by external body forces. Secondly, note that is the amount of dissipated energy in Ω in some time instant by a non-Newtonian fluid obeying the constitutive relation T = −pI + 1 0 S ε (λDv) dλ in Ω and the boundary condition (Tn) τ = − 1 0 s ε (λv) τ dλ on Γ N of the Navier-slip type, where T is the Cauchy stress tensor and p is the pressure. Finally, the term ε| • v| 2 represents the magnitude of certain inertial forces scaled by the parameter ε > 0.
We remark that ε has units of time and it models a certain future-time horizon that is still taken into account to get information about the present state. The causality is thus recovered only in the limit ε → 0 + ; we refer to [56] for a detailed discussion. The parameter ε has also a secondary role as it introduces additional dissipation via S ε and s ε , thus stabilizing the functional I ε . Similarly as in [56], we choose to work on the infinite time interval (0, ∞), which avoids prescription of the terminal condition for velocity. This would be otherwise necessary since the Euler-Lagrange equation corresponding to I ε is of the second order in time. To summarize, compared to the WIDE functional considered in [56], we make four important changes: 1) We remove the constraint v = 0 on ∂Ω and include friction on Γ N .
2) We modify the inertial term by considering This identity is unaffected by the boundary conditions and it is crucial for obtaining an energy estimate for the solution of Euler-Lagrange equation of I ε . It seems that our method works for the standard form of the inertial term | • v| 2 only if r > 3, leading to simpler outflow boundary conditions, without the corrector 1 2 |v| 2 n. 3) We develop the whole theory for certain non-Newtonian fluids with nonlinear dependence of the stress on strain, both in the bulk and on the boundary.

4)
We use a different stabilization term, which admits a certain physical interpretation. Note that the aim of stabilizing I ε is that of possibly obtaining information on |∂ t v| 2 and | rot v × v| 2 , starting from bounds on | Relation (1.10) just prescribes a condition of the Navier-slip type on the impermeable part of the boundary Γ N . The standard Navier-slip corresponds to the case where r = 2 and s is linear.
It is interesting to note that the first identity in (1.11) is strikingly similar to the stationary Navier-Stokes equation This is actually quite intuitive since the outlet is just an abstract boundary, where nothing should happen to the flow, and therefore the equation that holds there should be just a restriction of the equation that is satisfied by the flow in the bulk. Moreover, since (1.11) is obviously nonlinear, it cannot be easily treated as a constraint, unlike the usual do-nothing boundary conditions. This suggests that the outflow boundary conditions should be perceived as a special kind of PDEs on the outlet boundary. Then, since we work with weak solutions, it is not at all surprising that we are able to identify the outflow boundary conditions only in a weak sense. This has been also observed in [12] for a slightly different type of outflow boundary condition. The analogous remark actually applies also to the tangential part of the Navier-slip-type boundary condition (1.10).
Without the corrector 1 2 |v| 2 n, relation (1.11) is sometimes called constant traction boundary condition (cf. [38]), which is just the do-nothing boundary condition (2.1), but for the symmetric part of the velocity gradient. As we shall see, we have the freedom to replace Dv in the definition of I ε by other types of gradients, leading, for example, to boundary conditions involving ∇v rather than Dv. However, this is at expense of losing the physical meaning of the dissipation terms in I ε because of the possible failure of material frame indifference.
It has been observed experimentally in [29] that the correction 1 2 |v| 2 n on the outflow boundary has a positive effect on the flow characteristics. So far it has been unclear whether the boundary conditions such as (1.11) can be somehow derived, or if they are completely artificial, see the discussion in [34]. Our result suggests, that there actually might be a certain physical justification behind. On the other hand, in [58] it is argued that the outflow boundary condition which is just a version of (1.11) with the full velocity gradient (and renormalized pressure constant) is unphysical since it does not allow the Poiseuille flow in a straight cylindrical pipe (the streamlines are bent inwards near the outlet, see the figures in [58, Fig. 4 c)] or [70]). This is certainly true if the outlet is flat, however, one has the freedom of prescribing (1.11) on a curved outlet boundary (which is just an artificial interface, after all). This can partially compensate for the additional term 1 2 |v| 2 n, but not completely. In fact, choosing a suitable shape of the outlet boundaries should be probably seen as a part of the whole optimization problem. For simplicity however, in this work we assume that the outflow boundaries are fixed a apriori.

2.5.
On the pressure and the constants c i . Let us provide more insight into the redefinition of the pressure, that was made from (1.6) to (1.7). It is well known that, in case of no inflows/outflows, i.e., if v · n = 0 on ∂Ω, the pressure can be completely eliminated from the Navier-Stokes equations by the Leray projection. This happens because Ω ∇p · ϕ = 0 whenever div ϕ = 0 and ϕ · n = 0 on ∂Ω. Consequently, any substitution of the type p(t, x) =p(t, x, v(t, x)) does not change the problem. This is no longer true if v · n = 0, in general. Indeed, whenever the boundary condition involves pressure (which we know it will, cf. (1.11)), that pressure must correspond to the pressure appearing under the gradient operator in the Navier-Stokes equation.
It is a well known fact in the modeling of internal flows of incompressible fluids that the pressure is determined only up to a constant (in space). This is immediately seen from (1.8). We wish to point out that, although the outflow boundary condition (1.11) involves pressure, it is actually invariant with respect to the shift of p by a constant due to (1.12), and therefore the whole system (1.8), (1.9)-(1.11) retains the same property. This is in agreement with the physical intuition that varying the pressure by a same amount at all points of an incompressible fluid does not affect the flow, and there seems to be no reason why the presence of outflows should change this fact, since one can reasonably assume that the same fluid occupies the space also behind the outlet. Thus, one has always the freedom to impose a single additional condition on p, such as the value at a point, an integral average over some subdomain etc., provided that these quantities can be defined.
There is an interesting analogy between the pressure p and the constants c i . While p is a Lagrange multiplier to the constraint div v = 0, the constants c i are Lagrange multipliers corresponding to Γ i F v · n = F i , i = 0, . . . , n. Physically, constants c i represent certain generalized pressure drops (cf. [29]) and they can all be shifted by a common constant that is incorporated in the pressure p, without affecting the velocity v. It also seems possible to treat both p and c i as unknowns of the system and compute them implicitly by minimizing the functional over an enlarged function space without the constraints div w = 0 and Γ i F w·n = 0. The implicit methods based on J ε may turn out relevant in numerical implementations of the problem, since they give enhanced numerical stability and easier construction of the function spaces for the solution and test functions. However, this obviously leads to a saddle point problem and some further stabilization S is necessary to ensure even the existence of a minimum of J ε . Since it is hard to think of any physical justification behind S and there are many possible choices, we shall stick to the functional I ε , search for the solution in the spaces constrained by div v = 0, Γ i F v · n = F i , and then construct p and c i a posteriori from v.
In the remaining part of the paper, we provide a rigorous counterpart of the procedure outlined in the introduction.

Technical assumptions & definitions
In this section, we state precisely the hypotheses needed to prove our main results. The definition of required function spaces is given here as well.
3.1. Constitutive assumptions for S and s. Relations (2.3) and (2.5) are specified via assumptions on S and s, which then allow us to apply our results to a wide class of non-Newtonian fluids.
For the function S, we suppose that for all A, B ∈ R 3×3 sym and some r > 1. It is easy to see that under these conditions, the function S ε , defined in (2.2), fulfills for all A, B ∈ R 3×3 sym and some q > 1. Note that thanks to (3.2), there holds ∂ ∂A Similarly, we require that s satisfies the properties: (3.16) and, consequently, the function s ε defined in (2.4) fulfills and there holds ∂ ∂u for all u ∈ R 3 . The assumptions (3.7) and (3.18) could be omitted, however then, the identities (3.11) and (3.22), characterizing the dissipation potentials, become more complicated, see [22,Corollary.]. The above assumptions could be further generalized in many ways (e.g., by allowing anisotropic constitutive relations). We do not aim at maximum generality here and stick with the simple setting. The parameters r and q retain the same meaning throughout the whole paper and they will be eventually required to satisfy certain bounds. For a future use, let us collect the above assumptions into the hypothesis (H S ): The functions S, S ε , s, s ε satisfy ( Of course, this also leads to very poor information about the pressure in the Navier-Stokes equations (since one cannot even apply the L p -regularity for elliptic systems). Since our main result is of qualitative nature (identification of boundary conditions), this lack of information has limited effect in the analysis below.
Next, we observe that, if the flux rates through the boundaries Γ D , Γ 1 F and Γ 2 F are prescribed, one can use the incompressibility constraint div v = 0 in order to infer that also the net flux rate through Γ 0 F is determined. Indeed we have the formula Hence, we assume throughout that the sets Ω, Γ D , Γ N , Γ i F , i = 0, . . . , n, and Γ F are chosen in such a way that Moreover, we let Consequently, the connected components of the set Γ F are separated by (subsets of) Γ D ∪ Γ N . In particular, if |Γ D | = |Γ N | = 0, then ∂Ω = Γ F = Γ 0 F . Note that one could, in principle, allow existence of some neighboring outlets, but this would be physically counter-intuitive and present an unnecessary complication in the following.

Function spaces.
We classically denote the Lebesgue and Sobolev spaces by (L p (Ω; R d ), · p ) and (W 1,p (Ω; R d ), · 1,p ), 1 ≤ p ≤ ∞, d ∈ N, respectively. Next, let us define and then, for any 1 < p < ∞, we define p ′ := p p−1 and The spaces V 1,p and V 1,p div are equipped with the norm Note that any w ∈ V 1,p satisfies ∂Ω w · n = 0 (and its boundary values can be thus extended to a divergence-free vector field). In what follows, we often rely on the Korn-Poincaré inequality on V 1,p in the form  (3.23). Note that due to (3.23) and the trace theorem, the expression Dw p + w p;ΓN is an equivalent norm on V 1,p .
We shall also need the following special spaces of Lions-Magenes type defined on a an open connected subset G of ∂Ω. First, let us denote and then, we put The boundary conditions will be identified in the corresponding dual spaces.

Nonhomogeneous data.
In the problem under consideration, we wish to prescribe a nonhomogeneous initial datum u 0 : Ω → R 3 and boundary data v D : (0, ∞) × Γ D → R 3 , F i : (0, ∞) → R, i = 0, . . . , n. However, due to the nonlinearity of the problem, it seems difficult to formulate some explicit necessary conditions on this data that are needed for the existence of the corresponding global-in-time weak solution of Navier-Stokes equations, especially if the data are allowed to be timedependent. We refer to [7], [24], [25,IX.4], or [59] where this topic is (partially) treated using different approaches. Basically, one is asking whether it is possible to extend the boundary data v D , F i to a divergence free field v 0 of such regularity that the products (∇v)v 0 and (∇v 0 )v are under control. Obviously, this depends on many factors and there seems to be no agreement on how this extension should be constructed. To avoid this difficult question and, at the same time, to make our results applicable in real scenarios with (large) inflow and outflows, we assume that the data u 0 , v D , F i and f are admissible in the sense that for any ε > 0 there exists a function v ε 0 : (0, ∞) × Ω → R 3 with the following properties: and eventually also Admissible data indeed exist, at least in some simple scenarios. An obvious case is when v D and F i are time independent and u 0 ∈ V 1,q div . Then, we can simply choose v ε 0 (t) := u 0 for all t > 0 and ε > 0. It is easy to see that if v ε 0 exist, then there also exists a function v 0 ∈ X 2 div ∩X r div ∩ L ∞ (0, ∞; L 2 (Ω; R 3 )) with ∂ t v 0 ∈ (X 2 div ∩ X r div ) ′ and satisfying the same constraints as v ε 0 , to which a subsequence of {v ε } ε>0 converges weakly in the corresponding spaces (and there also holds v(0) = u 0 , div u 0 = 0).

Main result: variational resolution of Navier-Stokes equations with outflow boundary conditions
We are now in the position of stating the main result of the paper, which makes the oulined variational approach rigorous. and (H f ) be fulfilled. Then, for every ε > 0, the functional I ε attains a minimum v ε in the set U ε : If, in addition, the conditions (H B 0 ) and which is a limit of a not relabeled subsequence of {v ε } ε>0 in the sense that and which solves for all ϕ ∈ X 2 div ∩ X r div . Furthermore, let D ∈ L ∞ loc (0, ∞; R). Then, there exists a function and such that for all ψ ∈ C ∞ c ((0, ∞); V 1,r ).
Theorem 1 is proved in Section 5 below. We devote the rest of this section to several comments instead: (i) Since the setting of Theorem 1 enables to test (4.10) by v − v 0 , it is clear that we are dealing with a weak solution of the Leray-Hopf type. The total energy of the fluid can however be increased by the non-homogeneous and possibly unsteady inflow, therefore the standard form of the energy (in)equality is unavailable. This can be also seen in the proof of the uniform estimate below, where the physical energy of the (approximated) solution is estimated by the size of the data.
(iii) Note that if the test function ϕ from (4.10) vanishes on the whole ∂Ω, then which is the standard form of the convective term in the Navier-Stokes equations.
If ϕ does not vanish on ∂Ω, the term −(∇v) T v · ϕ contributes both to a pressure and to a boundary term as can be seen in (4.12).
(iv) In assumption (4.2), the parameters q, σ 4 , and ρ 4 correspond just to the stabilization of I ε and their values have no impact on the properties of (v, P ).
(v) The upper bounds r < 4 and q ≤ 3r ′ are easily removable. The inequality q ≤ 3r ′ just simplifies some dual spaces that are needed below. The case r ≥ 4 is sub-critical and the analysis becomes simpler as the stabilization of I ε by the fourth order dissipation terms is unnecessary, i.e., one can set σ 4 = ρ 4 = 0. Furthermore, if r > 4, we put σ 4 = ρ 4 = σ q = ρ q = 0 and the parameter q is omitted completely.
(vi) The lower bound r ≥ 11 5 is probably not optimal. However, in the case r ∈ ( 6 5 , 11 5 )\{2}, the identification of the weak limit of S(Dv ε ) becomes difficult. It is not obvious whether the standard methods (such as the L ∞ -or W 1,∞ -truncations) can be applied directly for our ε-approximation scheme, which is rather complicated.
Also, in the case r < 2, there is another difficulty of extending the boundary data without some kind of smallness condition on v D · n and F i , cf. e.g. [32] and [37].
(vii) Theorem 1 is valid, of course, for the classical Navier-Stokes system, where r = 2 and S, s are linear functions. However, the spaces for ∂ t v and for ϕ in (4.10) have to be modified in a standard way (it is no longer possible to consider v − v 0 as a test function in (4.10)). Also, it is apparent in the proof below (cf. (5.29)) that if r = 2, one has to replace the assumption v ε 0 L ∞ (0,∞;L ∞ ) ≤ C (recall (H B 0 )) by the requirement that the number ess sup can be made smaller than any given δ > 0. That v ε 0 can indeed be constructed in such a way follows from [25,(IX.4.43)] and it is closely related to the so-called extension condition of Leray and Hopf.
(viii) With usual modifications, the above theorem can be proved also in the d-dimensional, d ≥ 2, setting.
(iv) Starting with the definition of I ε , it is possible to replace everywhere the symmetric velocity gradient Dv by the full velocity gradient ∇v. This changes very little in the analysis below. Then, the analogue of Theorem 1 yields the outflow boundary condition −pn − 1 2 |v| 2 n + S(∇v)n = c i n on Γ i F (again, if the solution is sufficiently smooth). This is a nonlinear version of the donothing boundary condition with the dynamic pressure correction. One could even replace Dv by ∇v + a(∇v) T for any a = −1 and again, since div(∇v) T = ∇ div v = 0, this would not change anything except the resulting boundary conditions. (x) If r > 3, then (∇v)v · v becomes easily controllable and hence, we do not need to use (2.6). Thus, by replacing Making δ > 0 very small, this is as close as we can get to the classical do-nothing boundary condition (2.1), while retaining the global-in-time and large-data existence of a weak solution of the Navier-Stokes equations corresponding to S.
(xi) By subtracting the term ∞ 0 e − t ε ΓF g · v in the definition of I ε , we can include also a forcing through the outflow boundaries, that may arise as a reaction force to the fluid flow outside Ω. Then, if g ∈ L 2 (0, ∞; L 2 (Γ F ; R 3 )), this additional term can be handled as the one with f , thanks to the fact that we are able to control the whole Sobolev norms of the solution via (3.23). (In fact, one could be slightly more general here and observe that f , g ∈ L 2 L 2 + L r L r is sufficient to obtain Theorem 1.) This way, it is possible to derive a nonhomogeneous version of the boundary condition (4.16). On the other hand, it seems hard to imagine any kind of physical device that would yield exactly some given forcing g, apart from the case where g arises as a boundary value of some potential, in which case it can be included in the pressure. Thus, to ensure a viable physical interpretation, one should probably stick with the choice g = 0, which is what we do in this work.
(xii) Adding the term ∞ 0 e − t ε ΓN ε 2 |∂ t v| 2 in the definition of I ε leads (at least formally) to the so-called dynamic slip boundary conditions, which are studied in [28] for certain polymer melts. Applying the similar idea also on the outflow boundary Γ F , one can recover the conditions involving ∂ t v, that are studied analytically in [11] and which are quite popular in numerics, see, e.g., [40] and references therein. For simplicity of presentation (it would require major changes in the definitions of the function spaces below), we do not discuss these extensions in detail here.
(xiii) The distribution p := ∂ t P from the above theorem is called pressure. Note that the boundary conditions are identified only in an averaged sense over time. This is to avoid the fact that div T, where T := −pI + S(Dv), is not a function, in general, and hence the trace of Tn cannot be defined consistently. Actually, this issue arises already for the non-stationary Stokes system with, e.g., the Navier-slip boundary condition if the domain or the initial data are rough and no additional regularity of solution, besides the energy estimate, is available. On the other hand, one can read from (4.12) that div T ψ is an integrable function, hereby providing a possibility to define the object T ψ n in a way which is compatible with the common meaning of a trace (a limit of restrictions of smooth functions) and with the corresponding integration-by-parts formula, that is essential for the identification of (4.15) and (4.16). For a different resolution of this issue, see [12]. It is clear that if T, div T ∈ L 1 (Q ∞ ; R 3 ), then (·) ψ averaging is unnecessary and conditions (4.15), (4.16) hold a.e. in (0, ∞). The functions ξ i always exist (see the end of the proof of Theorem 1 below) and their only purpose is to get around the fact that the functional T ψ n cannot be applied to n, for a general Lipschitz domain.

Proof of the main result
The proof of Theorem 1 is developed throughout the section and and is divided in subsequent sections and structured on a series of lemmas. In what follows, the number ω ∈ (0, 1 2 ) is systematically used as an auxiliary parameter (arising, e.g., from the Young inequality), which is eventually chosen sufficiently small.

Existence of minima.
Let us consider the problem of minimizing I ε over the space U ε . The following lemma tells us that this problem is solvable. , and (H f ) are satisfied. Then, the functional I ε attains a minimum v ε in U ε . Moreover, there exists a constant C > 0, depending only on the data and Ω, such that Proof. Let us start by deriving some preliminary estimates that are used below for estimation of certain terms of I ε .
Let {v k } ∞ k=1 ⊂ U ε be a minimizing sequence for I ε , satisfying also I ε (v k ) < C for all k ∈ N. By using (5.5), (5.6), (5.9), (5.10) and then choosing ω > 0 sufficiently small, we get Using this estimate on (0, T ), T > 0, we can apply the Aubin-Lions lemma with the compact embedding W 1,2 (Ω; R 3 ) ֒→ L 2 (Ω; R 3 ) to deduce that a nonrelabeled subsequence of v k converges to v ε ∈ U ε strongly in L 2 (0, T ; L 2 (Ω; R 3 )) and ∇v k converges weakly to ∇v ε . Therefore, the product rot v k × v k , which is bounded in L 2 (0, T ; L 2 (Ω; R 3 )) (recall (5.4)), converges weakly (up to another subsequence) to the correct limit rot v ε × v ε . Hence, by the weak lower semi-continuity of norms, we obtain and, taking the limit T → ∞, we see that v ε is a minimum of I ε .

5.2.
Euler-Lagrange equations. Next, we derive the Euler-Lagrange equation of I ε , which we state in three different ways: (5.13) is useful for deriving most of the ε-uniform estimates, (4.1) is convenient for the weak limit identification, while (5.15) is important to get an ε-uniform estimate of ∂ t v ε .
Proof. Let ϕ ∈ V 0 . Since V 0 ⊂ U 0 , we have, for all s ∈ R, that v ε + sϕ ∈ U ε and it is easy to see that also I(v ε + sϕ) < ∞. Let us define By using the fundamental theorem of calculus in the form g(s) − g(0) = s An analogous identity holds, of course, for the function s ε . Hence, also using we find that Dividing this by s > 0, using I ε (ϕ) < ∞ and the continuity of S ε and s ε to take the limit s → 0 + , and then doing the analogous procedure for s < 0 leads to (5.13).
If ψ ∈ V c , then also ϕ := e To prove (5.15), note first using Hölder's inequality and the fact that v ε ∈ U ε that the functionals N 1 and Av ε are well-defined on V 1, 2q q−2 , V 1, q q−3 , V 1,4 , V 1,q and V 1,r , respectively, for almost every time, where the space V 1,q is the smallest one of these. Indeed, this a consequence of the inequality that follows (5.14). Hence, if we rewrite (4.1) as ), leading to (5.15).

5.3.
A priori estimates. Now we examine the limit of the sequence of minimizers constructed in the previous section. First, we need to derive an ε-uniform estimate for v ε , which is the key technical point of the paper.
Proof. To get an estimate on (0, ε), we notice that e − t ε ≥ e −1 on that interval, hence we can just apply (5.1), use Young's inequality and (3.23) to eliminate f and then estimate the exponential from below. This leads to Next, we let T > 0 and choose ϕ : Note that then ϕ ∈ V 0 and ϕ(0) = 0. Hence, we obtain where we used (2.6) in order to eliminate the term (∇v ε )v ε · v ε , which would otherwise be impossible to control for general boundary conditions and r ≤ 3. Next, we apply the identity Our aim is to estimate the terms on the right-hand side by the ω-Young inequality, so that the part containing v ε is sufficiently small and the other (large) part depends only on v ε 0 , for which we can use assumptions (H A 0 ), (H B 0 ), (H f ) and include these terms into a constant C (or C(ω)). Let us proceed term by term and observe first, using (3.23), that The products ε| rot v ε × v ε | 2 are already estimated in (5.4). Note that the assumptions on σ 4 and ρ 4 are needed to absorb the term 1 4 ε| rot v ε × v ε | 4 which is the only term on the right-hand side of (5.24) that does not involve data. Next, for the term with S ε , we use (3.4) and Young's inequality to get and, similarly, we also obtain To handle the term containing f , we proceed analogously as in (5.6), leading to Regarding the last two terms in (5.24), we use (3.23) to get and then the Bochner version of integration by parts formula to write where in the last inequality we also used (5.21) to estimate the integral over (0, ε) by a constant. If we apply the inequality , t > 0, on the left-hand side of (5.30), we see that the integral on the right-hand side of (5.30) gets absorbed for ω sufficiently small, leading to Putting this information together with (5.21) and then taking the essential supremum over T > 0, we arrive at and also at By virtue of (3.23), estimate (5.31) is equivalent to (5.17). Moreover, the information ε 2) which, together with (5.32), yields (5.18) through the Young inequality. Moreover, by a similar estimate to (5.29), we also immediately obtain (5.19) from (5.17).
A direct calculation verifies that, in the considered range of r, we have 1 < z ≤ 33 4 ≤ 3r 3−r , hence the embeddings hold true, from which we deduce This together with (5.34) gives for any r ≥ 11 5 . Next, we estimate • v ε by applying a similar method as in [56]. Estimates (5.17)-(5.19) proved thus far and the Hölder inequality show that the functionals N 1 ε v ε and Av ε , defined in Lemma 3, are bounded in the following sense: Then, we use (5.15) to express • v ε as a temporal convolution with the kernel K(t) := ε −1 e t ε χ t≤0 . This leads to which is understood as an identity in the space L 1 loc (0, ∞; V −1,q ′ div ). Using the properties of convolution, (5.36)-(5.40), and recalling (5.14), the right-hand side of (5.41) is a continuous linear functional in the space X 2 ∩ X q . Identity (5.41) thus gives 5.42) and in combination with (5.35) and the embedding X 2 div ∩ X q div ֒→ X r , this proves (5.20).
It is clearly seen in the proof above that the assumption r ≥ 11 5 is used only to show that the convective term is a bounded functional on X r , uniformly with respect to ε > 0. This information is useful later when taking the limit ε → 0 + . Otherwise, it is important that r ≥ 2, because then one can absorb the term (5.29) and remark (vi) above. It is apparent that the term ε|Dv ε | 4 is needed to control terms related to the convective term. The role of the higher order stabilization ε|Dv ε | q , q > 4, is later clarified while taking the limit ε → 0 + .
Comparing with the usual existence theories for Navier-Stokes equations, one may wonder why we need a super-linear growth also in the boundary terms on Γ N . We recall that we do not want to impose additional geometrical assumptions on Ω. In this case, there will always be a term on the right-hand side of (3.23) that controls the overall speed of the flow (to have just w 1,p ≤ c ′ p Dw p , one would need to exclude domains that are "too special" such as axisymmetric domains, parallel plates etc.) As opposed to usual existence theories, we cannot choose this term to be ( Ω |w|) p since we do not know that v ε L ∞ (0,∞;L 2 ) is bounded a priori. In our case, this information needs to be carefully deduced from the Euler-Lagrange equations by testing with a solution (minimum), but this generates many terms without a sign, especially in the case with nonhomogeneous data, as can be seen in the proof above. It thus seems natural to take instead into consideration the fact that the fluid loses energy also due to friction on Γ N . But then (3.23) indirectly requires the scalings of S and s to be compatible, explaining the same nonlinear growth. The whole situation would simplify in the case we considered I ε for the full velocity gradient, since then one does not need the Korn inequality. Nevertheless, even for the Poincaré inequality to hold in the form w 1,p ≤ c ′′ p ∇w p , certain domains have to be ruled out.

Pressure reconstruction.
At this point, the most of the work leading to Theorem 1 is done. Before proceeding with its proof, let us state one more auxiliary result, that is used in the pressure construction. Let us define L p 0 (Ω; R) := {f ∈ L p (Ω; R) : Ω f = 0}. Up to the boundary conditions, the following proposition is very standard.
To see more explicitly that the boundary conditions encoded in V 1,p do not cause any difficulties, one can show that the auxiliary problem for all ϕ ∈ V 1,p admits an unique solution (u, q) ∈ V 1,p × L p ′ 0 (Ω; R) whenever g ∈ V −1,p ′ , which then obviously gives (5.44) if (5.43) holds. This is nothing but the weak formulation of the nonlinear problem div u = 0, − div(|∇u| p−2 ∇u) + |u| p−2 u + ∇q = g in Ω, where c i are the Lagrange multipliers to the constraints Γ i This condition is again an immediate consequence of W 1,p 0 (Ω; R) ⊂ V 1,p , the fact that the norms · 1,p and ∇· p are equivalent on W 1,p 0 (Ω; R 3 ) and the standard inf-sup condition for the pair (W 1,p 0 (Ω; R 3 ), L p ′ 0 (R)) (to be found in various forms in the works by O. A. Ladyzhenskaya, J.-L. Lions, E. Magenes, I. Babuška, J. Nečas, or F. Brezzi), which can be proved by applying the Bogovskii operator to |q| p ′ −2 q − 1 |Ω| Ω |q| p ′ −2 q ∈ L p 0 (Ω; R).
Note that Proposition 5 works for a general Lipschitz domain Ω (actually only the local cone property is needed), which is desirable in our application, cf. Figure 1. This contrasts with other methods of constructing q, such as the Helmholtz decomposition or the L p -theory for the Stokes system that require some regularity of ∂Ω, cf. [8] and references therein. 5.5. Passage to the limit as ε → 0 + . Let {v ε } ε>0 be a sequence of minimizers to I ε , which exist due to Lemmas 2 and 3.
By the properties of the trace operator (see [16,Corollary 1.13.]), it is standard to show that the trace of v ε actually converges strongly to the trace of v on Γ N (proving the second part of (4.8)) and then, by the continuity of s, this necessarily means that z = s(v).
To prove (4.10), it remains to identify the weak limit S. To this end, we take advantage of the fact that in the considered case r ≥ 11 5 , the function v, after a correction of boundary values, is an admissible test function in (5.56). Let v δ be a fixed element of the approximating sequence {v ε } ε>0 and let 0 ≤ η ∈ C ∞ c ((0, T )), T > 0. Next, we observe that and are well-defined and finite quantities. Moreover, using the convergence results (4.6)-(4.9), (5.47)-(5.54), Lemma 4 and the property v Further, using ϕ := (v − v δ )η ∈ X 2 div ∩ X r div as a test function in (5.56) leads to Next, we use ψ := (v ε − v δ )η ∈ X 2 div ∩ X q div in (4.1), (5.23), Young's and Hölder's inequalities, q > 4, and (5.2), giving Now, we take the limes superior of this inequality and on the left-hand side we use that v ε → v strongly in L 2 (Q T ; R 3 ) (by interpolation and Vitali's theorem), the inequality S ε (Dv ε ) · Dv ε ≥ S(Dv ε ) · Dv ε , and in the boundary term we use (4.8) and Fatou's lemma. This way, we get Comparing this with (5.60) immediately leads to lim sup Hence, by the monotonicity of S, we get, for any W ∈ L r (Q T ; R 3×3 sym ), that Choosing now W = Dv + λU , U ∈ L r (Q T ; R 3×3 sym ) and dividing by λ > 0 yields 0 ≤ QT (S(Dv + λU ) − S) · U η and, consequently, using the continuity of S to take the limit λ → 0 + , we arrive at 0 ≤ QT (S(Dv) − S)η · U.
Since U is arbitrary, we deduce that S(Dv)η = Sη a.e. in Q T , but since η and T are also arbitrary, we conclude that S(Dv) = S a.e. in Q ∞ and (4.10) is proved.
In the next step, we prove (4.12) by constructing a pressure in (4.10). Since the test functions from (4.10) must vanish on Γ D , we follow the same construction of pressure as in [71, Theorem 2.6.], but only partially, since we do not need a pressure decomposition here.
Recalling (5.66) and (3.25), this remains valid for all w ∈ W 1 r ′ ,r 0,n (Γ N ; R 3 ), proving (4.15). Next, let 0 ≤ i ≤ n and choose w ∈ W 1,∞ 0 (Γ i F ; R 3 ) such that Γ i F w · n = 0, extend this function by zero to whole ∂Ω, and use (5.68) to deduce where we abbreviated The restriction of w on Γ i F is a consequence of the prescribed net fluxes and of the incompressibility of v, recall the definition of V 1,r . Since Γ i F is locally a graph of a Lipschitz function, there exists a vector field ξ i ∈ W 1,∞ (Γ i F ; R 3 ) such that ξ i · n > 0 on Γ i F . Moreover, it is clear that ξ i can be chosen in a way that ξ i ∈ W 1,∞ 0 (Γ i F ; R 3 ) and Γ i F ξ i · n = 1. Let z ∈ W 1,∞ 0 (Γ i F ; R 3 ). Then, the function can be used in (5.69), leading to Since z ∈ W 1,∞ 0 (Γ i F ; R 3 ) was arbitrary and n ∈ L ∞ (Γ i F ; R 3 ), we deduce (4.16) by virtue of (3.24) and (5.66).
The bounds r ≥ 11 5 and q > 4 are evidently used only to identify that the weak limit of S ε (Dv ε ) is S(Dv). The assumption r ≥ 11 5 greatly simplifies the identification procedure since then the function v can be used (after minor corrections) in (4.10) as a test function. But it is unlikely that this bound for r is necessary since in the mathematical theory of non-Newtonian fluids, more refined arguments are known for the limit identification. Unfortunately, the application of the methods of either [71], or [14] to our ε-approximation scheme seems not straightforward. Therefore, the case 6 5 < r < 11 5 , r = 2 is left open. However, this drawback seems not so significant in the view of the fact that our method works for p-fluids with p sufficiently large that can approximate a flow of an r-fluid for arbitrary r > 6 5 , see the proof of [14, Theorem 3.1.], effectively avoiding the ε-approximation.