<figure><img src="image_1.png" /></figure>

RESEARCH ARTICLE

<figure><img src="image_2.png" /></figure>

# Stability of a continuous/discrete sensitivity model for the Navier-Stokes equations

N. Nouaime<fnref n="1" />,<fnref n="2" /> ID

B. Després<fnref n="2" />

M. A. Puscas<fnref n="1" />

C. Fiorini<fnref n="3" />

**Correspondence**

N. Nouaime, Université Paris-Saclay,
CEA, Service de Thermo-hydraulique et
de Mécanique des Fluides, 91191
Gif-sur-Yvette, France.

Email: nathalie.nouaime@cea.fr

**Abstract**

This work presents a comprehensive framework for the sensitivity analysis of the Navier-Stokes equations, with an emphasis on the stability estimate of the discretized first-order sensitivity of the Navier-Stokes equations. The first-order sensitivity of the Navier-Stokes equations is defined using the polynomial chaos method, and a finite element-volume numerical scheme for the Navier-Stokes equations is suggested. This numerical method is integrated into the open-source industrial code TrioCFD developed by the CEA. The finite element-volume discretization is extended to the first-order sensitivity Navier-Stokes equations, and the most significant and original point is the discretization of the nonlinear term. A stability estimate for continuous and discrete Navier-Stokes equations is established. Finally, numerical tests are presented to evaluate the polynomial chaos method and to compare it to the Monte Carlo and Taylor expansion methods.

**KEYWORDS**

finite element-volume scheme, Navier-Stokes equations, polynomial chaos method

# 1 | INTRODUCTION

The use of numerical simulation has grown significantly in various fields, including nuclear engineering, where it has become an essential tool for ensuring safety and analyzing accidents. In this context, a modern approach is to simulate systems of Partial Differential Equations (PDEs), such as the Navier-Stokes equations, using the Polynomial Chaos Methodology (PCM) to account for uncertainties that are naturally present in physical models. With PCM, one is able to compute both the main model's solution and additional sensitivity variables. In this article, we apply the PCM to the Navier-Stokes equations, obtaining a system of equations that describes the dynamics of the physical problem and its sensitivity to changes in parameters. This system is twice the size of the original system and fully coupled. In this work, we propose a numerical scheme for this system, which is implemented in the open-source industrial platform TrioCFD, and we analyze the stability properties in quadratic norm, both for its continuous and discrete version. We also introduce a change of variable to partially decouple the obtained system. This allows for a less intrusive and faster to implement method.

The PCM is one possible approach to tackle Uncertainty Quantification (UQ) problems. PCM has been used to solve numerous problems in various domains. It is focused on propagating uncertainties from inputs to outputs of a numerical model and provides an efficient and accurate way of dealing with uncertainty.<fnref n="1" />⁻⁵ This task can be performed in many different ways, and it may vary depending on factors such as the nature of the model, the amplitude of perturbations, and

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2024 The Author(s). *International Journal for Numerical Methods in Fluids* published by John Wiley & Sons Ltd.

<footnote n="1">¹Université Paris-Saclay, CEA, Service de Thermo-hydraulique et de Mécanique des Fluides, Gif-sur-Yvette, France</footnote>

<footnote n="2">²Sorbonne Université, CNRS Université de Paris Laboratoire Jacques-Louis Lions (LJLL), Paris, France</footnote>

<footnote n="3">³M2N, Conservatoire National des Arts et Métiers, Paris, France</footnote>
whether they are deterministic or stochastic. In this work, we only consider systems governed by PDEs and we estimate the uncertainty of the model's solutions caused by uncertain input parameters. PCM allows the treatment of stochastic variables that can be described by their Probability Density Functions (PDF). More precisely, these variables are expressed as a linear combination of orthogonal polynomial functions of normalized random variables.

Norbert Wiener introduced the Polynomial Chaos Expansion (PCE) in 1938, in Reference 6, employing Hermite polynomials to model stochastic processes involving Gaussian random variables. This technique can be considered the original PCE method. Xiu and Karniadakis 7,8 then generalized this method to include other polynomials. Each polynomial has an optimal random distribution attached to it, giving a fast convergence rate. This article applies the PCM to the incompressible Navier-Stokes equations with normally distributed uncertain parameters. The Navier-Stokes equations are a system of PDEs which describe the motion of incompressible fluids. They consist of two equations, the conservation of mass equation, and the conservation of momentum equation. The equations are nonlinear and involve second-order derivatives of the fluid velocity, making them difficult to solve analytically. Nevertheless, they are widely used in engineering and science to model the behavior of fluids in various applications, from predicting the weather to designing aircrafts. 9,10 The study of the Navier-Stokes equations is a central topic in fluid dynamics, and remains an active area of research today. In a previous work, 11 the Navier-Stokes equations were considered, and the first-order sensitivity of the Navier-Stokes equations was obtained using the Taylor expansion. 12,13 In this approach, the sensitivity is defined as the derivative of the model's output with respect to the parameters of interest. 14 It can be used, among other things, to estimate of the variance of the solution of the Navier-Stokes equations when there are uncertain parameters. Confidence intervals can be provided, based on the estimated variance.

In References 11 and 15, the authors also discretize the sensitivity of the Navier-Stokes equations using the Finite Element-Volume (FEV) scheme 16–18 and in Reference 11 a stability estimate for the continuous sensitivity of the Navier-Stokes equations is provided. However, the stability estimate of the discrete Navier-Stokes and the discrete sensitivity equations were not established. The FEV method was originally formulated by Emonot 16 with applications to the Laplacian and Stokes problems. This method, called the FEV $P_{NC}^1 \setminus P^0$, is a modification of the Crouziex-Raviart element. 19 Heib, in Reference 18, showed that the element proposed by Emonot had several defects. Based on this observation, he proposed a new discretization by modifying the $P_{NC}^1 \setminus P^1 + \text{Bubble}$ element, for which, however, there are no convergence results for the pressure. Finally, Heib defined and applied the $P_{NC}^1 \setminus P^0 + P^1$ element and applied it to Stokes problems in Reference 18. For the $P_{NC}^1 \setminus P^0 + P^1$ element, there is a convergence result 18 for the velocity and the pressure of the Stokes equations. This numerical scheme is integrated in the open-source fluid dynamics simulation software, TrioCFD, 20,21 developed by the CEA.

The main original contributions of this work are the following:

1. We introduce a change of variable, thanks to which the Navier-Stokes equations and their sensitivity equations are partially decoupled Equation (20)–(21) (i.e., the sensitivity system's solution depends on the state system's one, but not the other way around). This new formulation is a priori simpler to use for discretization and less intrusive. We show the quadratic stability of the new system in Proposition 2.

2. We present the FEV scheme that we will use to discretize the sensitivity system. Preparatory to further development for the analysis of the sensitivity system, we show a new stability estimate for the discrete Navier-Stokes equations in Appendix A.

3. Then, we use all the previous results to obtain the main theoretical result of this work, which is the stability estimate for the FEV scheme for the first-order sensitivity of the Navier-Stokes equations, in Theorem 1.

4. Finally, we present some numerical results when we have an uncertain boundary condition. We compute the mean, the variance, and the confidence intervals for the velocity and the pressure.

This article is organized as follows. In Section 2, the Navier-Stokes equations with homogeneous Dirichlet boundary conditions are presented. Then, the polynomial chaos method is presented, and the derivation of the first-order sensitivity equations, according to this method, is detailed in Section 3. Section 4 provides a detailed explanation of the FEV method, especially the $P_{NC}^1 \setminus P^0 + P^1$ discretisation. In Sections 5, the Navier-Stokes equations are discretized according to the FEV method, and a useful and important property of the discrete trilinear term is proved, then a stability estimate is established for the discrete system. We apply in Section 6 the FEV method to the sensitivity of the Navier-Stokes equations and we prove that this discretization is stable in $L^2(0, T; H^1(\Omega)) \cap L^\infty(0, T; L^2(\Omega))$. To illustrate the method presented in the previous sections, we will show some numerical results in Section 7 and compare them with those obtained by other
sensitivity analysis methods such as Taylor expansion and Monte Carlo. Finally, in appendices A and B, the stability estimate of the discrete Navier-Stokes equations and the continuity of the linear form in the context of the FEV method are proved respectively.

## 2 | NAVIER-STOKES EQUATIONS

Let $\Omega \in \mathbb{R}^2$ be a domain that has a Lipschitz-continuous boundary, $\mathbf{X} = H_0^1(\Omega)^2 = H_0^1(\Omega) = \{\mathbf{v} \in H^1(\Omega) : \mathbf{v}|_{\partial\Omega} = \mathbf{0}\}$ and $W = L_0^2(\Omega) = \{p \in L^2(\Omega) : \int_{\Omega} p \, dx = 0\}$. The scalar product defined on $\mathbf{X}$ is $(\cdot, \cdot)_\mathbf{X}$ and $(\cdot, \cdot)$ is the scalar product defined on $L^2$. The incompressible Navier-Stokes system consists in finding $(\mathbf{u}, p)$ such that:

$$ \left\{ \begin{array}{ll} \partial_t \mathbf{u}(\mathbf{x}, t) - \nu \Delta \mathbf{u}(\mathbf{x}, t) + (\mathbf{u}(\mathbf{x}, t) \cdot \nabla) \mathbf{u}(\mathbf{x}, t) + \frac{1}{\rho} \nabla p(\mathbf{x}, t) = \mathbf{f}(\mathbf{x}), & \Omega, t > 0, \end{array} \right. \quad (1a) $$

$$ \nabla \cdot \mathbf{u}(\mathbf{x}, t) = \mathbf{0}, \qquad \Omega, t > 0, \quad (1b) $$

$$ \mathbf{u}(\mathbf{x}, 0) = \mathbf{0}, \qquad \Omega, t = 0, \quad (1c) $$

$$ \mathbf{u}(\mathbf{x}, t) = \mathbf{0} \qquad \text{on } \Gamma = \partial\Omega, t > 0. \quad (1d) $$

where $\mathbf{u} = (u^x, u^y)$ is the velocity, $p$ the pressure, $\mathbf{f}$ the external force, $\nu$ the kinetic viscosity and $\rho$ the density of the fluid. The first equation models the conservation of the momentum and the second one the conservation of the mass. For the sake of simplicity, only homogeneous Dirichlet boundary conditions are considered. A recent work in Reference 11 deals with physically motivated non-homogeneous Dirichlet boundary conditions, that we will use in the numerical section. The weak formulation of the Navier-Stokes equations (1) is written as follows

$$ \left\{ \begin{array}{ll} (\partial_t \mathbf{u}, \mathbf{v}) + a(\mathbf{u}, \mathbf{v}) + \frac{1}{\rho} b(\mathbf{v}, p) + t(\mathbf{u}, \mathbf{u}, \mathbf{v}) = l(\mathbf{v}), & \forall \mathbf{v} \in \mathbf{X}, \\ b(\mathbf{u}, q) = 0, & \forall q \in W, \end{array} \right. \quad (2) $$

with $a(\mathbf{u}, \mathbf{v})$ the bilinear form on $\mathbf{X} \times \mathbf{X}$

$$ a(\mathbf{u}, \mathbf{v}) = \int_{\Omega} \nabla \mathbf{u} : \nabla \mathbf{v} \, d\mathbf{x}, \quad (3) $$

$b(\mathbf{v}, p)$ the bilinear form on $\mathbf{X} \times W$

$$ b(\mathbf{v}, p) = \int_{\Omega} \nabla p \, v \, d\mathbf{x} = \int_{\Omega} \nabla \cdot \mathbf{v} p \, d\mathbf{x} + \int_{\partial\Omega} \mathbf{v} \cdot \mathbf{n} p \, dS = \int_{\Omega} \nabla \cdot \mathbf{v} p \, d\mathbf{x}, \quad (4) $$

$l(\mathbf{v})$ the linear form on $\mathbf{X}$

$$ l(\mathbf{v}) = \int_{\Omega} \mathbf{f} \cdot \mathbf{v} \, d\mathbf{x}, \quad (5) $$

and $t(\mathbf{u}, \mathbf{v}, \mathbf{w})$ the trilinear form on $\mathbf{X} \times \mathbf{X} \times \mathbf{X}$

$$ t(\mathbf{u}, \mathbf{v}, \mathbf{w}) = \int_{\Omega} [(\mathbf{u} \cdot \nabla)\mathbf{v}] \cdot \mathbf{w} \, d\mathbf{x}. \quad (6) $$

Throughout the article, we will use the following classical result.22,23

**Proposition 1.** For all $\mathbf{u}, \mathbf{v}, \mathbf{w} \in \mathbf{X}$ such that $\nabla \cdot \mathbf{u} = 0$, one has

$$ t(\mathbf{u}, \mathbf{v}, \mathbf{v}) = 0 \text{ and } t(\mathbf{u}, \mathbf{v}, \mathbf{w}) = -t(\mathbf{u}, \mathbf{w}, \mathbf{v}). $$
## 3 | POLYNOMIAL CHAOS METHOD

Sensitivity analysis studies how changes in a model's inputs affect its outputs. There are many methods to perform a sensitivity analysis, like Monte Carlo simulations,²⁴ regression analysis,²⁵ PCM,³ and the traditional approach is based on the Taylor series expansion.¹¹,¹⁴,²⁶ This article focuses on the PCM. There are two types of PC approaches: intrusive and non-intrusive. Our work is based on the Intrusive Polynomial Chaos (IPC) approach. This approach replaces all dependent variables in the Navier-Stokes equations by their PCE of order n to model the uncertainty. The resulting equations, once projected onto orthogonal polynomials, yield n + 1 times the number of deterministic equations. The obtained system can be solved using the same numerical methods as the original deterministic equations. This projection also provides a convenient representation of the variability in the model output with respect to the inputs, giving rise to a new model commonly referred to as the sensitivity model. The sensitivity model, in essence, allows us to understand how changes in the input variables affect the output, providing valuable insights into the system's behavior. This section introduces the PCM and presents the first order sensitivity of the Navier-Stokes equations, Equation (1), using this method.

### 3.1 | Polynomial chaos expansion and polynomial computation

Let $Y(\mathbf{x}, t; a)$ be a physical variable. Typically in our case, this physical variable can be the horizontal or vertical velocity or the pressure. The variable depends on $\mathbf{x}$, $t$ and is assumed to depend on an additional uncertain parameter $a$, normally distributed with $\mu$ the average and $\sigma$ the standard deviation, $a \sim N(\mu, \sigma^2)$. In PCM, the uncertain variables can be decomposed on the basis of complete orthogonal polynomials, the so-called PCE.¹,²,²⁷ The variable $Y(\mathbf{x}, t; a)$ can then be expressed by its PCE

$$Y(\mathbf{x}, t; a) = \sum_{i=0}^{n} Y_i(\mathbf{x}, t) \psi_i(a). \qquad (7)$$

The unknowns, $Y_i$, are deterministic coefficients and represent the random mode $i$ of the physical variable $Y$, and $\psi_i$ are the orthogonal polynomials of degree $i$. In the present work, only first-order sensitivity is considered that is, $n = 1$. Orthogonality means that

$$\langle \psi_i, \psi_j \rangle = \langle \psi_i, \psi_i \rangle \delta_{ij}, \qquad (8)$$

where $\langle \cdot, \cdot \rangle$ denotes the inner product

$$\langle \psi_i, \psi_j \rangle \equiv \int w(a) \psi_i(a) \psi_j(a) da, \qquad (9)$$

where $w$ is the weighting function and $\delta_{ij}$ is the Kronecker delta. If $\langle \psi_i, \psi_j \rangle = 1$, the polynomials are called orthonormal. The optimal polynomials $\psi_0(a)$ and $\psi_1(a)$ for a normally distributed random variable $a$ are given here below. For this purpose the weighting function $w(a)$ is considered to be the PDF of the Gaussian distribution $w(a) = \frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(a-\mu)^2}{2\sigma^2}}$ with $\int_R w(a)da = 1$.

According to the inner product Equation (9), the orthonormal polynomials are the following

$$\psi_0(a) = 1 \text{ and } \psi_1(a) = \frac{a - \mu}{\sigma}.$$

*Remark 1.* In general, Hermite, Legendre, Laguerre, Jacobi, and generalized Laguerre orthogonal polynomials are used for modeling the effect of uncertain variables described by standard normal distribution, uniform, exponential, beta, and gamma probability distributions, respectively.⁸

### 3.2 | Sensitivity of the Navier-Stokes equations

To obtain the first-order sensitivity of the Navier-Stokes system (1), first, the variables, that is, the velocity $\mathbf{u}$, the pressure $p$, and the external force $\mathbf{f}$ need to be expressed by their PCE:
$$ \mathbf{u}(\mathbf{x}, t; a) = \mathbf{u}_0(\mathbf{x}, t)\psi_0(a) + \mathbf{u}_1(\mathbf{x}, t)\psi_1(a) = \mathbf{u}_0(\mathbf{x}, t) + \left(\frac{a-\mu}{\sigma}\right)\mathbf{u}_1(\mathbf{x}, t), $$

$$ p(\mathbf{x}, t; a) = p_0(\mathbf{x}, t)\psi_0(a) + p_1(\mathbf{x}, t)\psi_1(a) = p_0(\mathbf{x}, t) + \left(\frac{a-\mu}{\sigma}\right)p_1(\mathbf{x}, t), $$

$$ \mathbf{f}(\mathbf{x}, t; a) = \mathbf{f}_0(\mathbf{x}, t)\psi_0(a) + \mathbf{f}_1(\mathbf{x}, t)\psi_1(a) = \mathbf{f}_0(\mathbf{x}, t) + \left(\frac{a-\mu}{\sigma}\right)\mathbf{f}_1(\mathbf{x}, t). $$

Second, these quantities are all inserted in Equation (1). Then each equation of this system is multiplied by $\psi_i$ (for $i = 0$ and $i = 1$) and to obtain the sensitivity equations the inner product is used.

For $i = 0$, the conservation of the momentum Equation (1a) becomes

$$ \int_{\Omega} \left[ \partial_t (\mathbf{u}_0(\mathbf{x}, t)\psi_0(a) + \mathbf{u}_1(\mathbf{x}, t)\psi_1(a)) - v\Delta(\mathbf{u}_0(\mathbf{x}, t)\psi_0(a) + \mathbf{u}_1(\mathbf{x}, t)\psi_1(a)) \right. \\ \left. + ((\mathbf{u}_0(\mathbf{x}, t)\psi_0(a) + \mathbf{u}_1(\mathbf{x}, t)\psi_1(a)) \cdot \nabla)(\mathbf{u}_0(\mathbf{x}, t)\psi_0(a) + \mathbf{u}_1(\mathbf{x}, t)\psi_1(a)) \right. \\ \left. + \frac{1}{\rho}\nabla(p_0(\mathbf{x}, t)\psi_0(a) + p_1(\mathbf{x}, t)\psi_1(a)) \right] \psi_0(a)w(a) da \\ = \int_{\Omega} (\mathbf{f}_0(\mathbf{x}, t)\psi_0(a) + \mathbf{f}_1(\mathbf{x}, t)\psi_1(a))\psi_0(a)w(a) da \quad \forall(\mathbf{x}, t). \tag{10} $$

The Equation (10) implies that for all $\mathbf{x} \in \Omega$ and $t \ge 0$,

$$ \partial_t \mathbf{u}_0(\mathbf{x}, t) - v \Delta \mathbf{u}_0(\mathbf{x}, t) + (\mathbf{u}_0(\mathbf{x}, t) \cdot \nabla) \mathbf{u}_0(\mathbf{x}, t) + (\mathbf{u}_1(\mathbf{x}, t) \cdot \nabla) \mathbf{u}_1(\mathbf{x}, t) + \frac{1}{\rho} \nabla p_0(\mathbf{x}, t) = \mathbf{f}_0. $$

The application of the PCM to the mass conservation Equation (1b) gives

$$ \int_{\Omega} [\nabla \cdot (\mathbf{u}_0(\mathbf{x}, t)\psi_0(a) + \mathbf{u}_1(\mathbf{x}, t)\psi_1(a))]\psi_0(a)w(a) = 0 \quad \forall(\mathbf{x}, t). \tag{11} $$

Therefore for all $\mathbf{x} \in \Omega$ and $t \ge 0$, $\nabla \cdot \mathbf{u}_0(\mathbf{x}, t) = 0$. The application of this method on the initial condition (1c) and the boundary condition (1d) respectively gives

$$ \int_{\Omega} [\mathbf{u}_0(\mathbf{x}, 0)\psi_0(a) + \mathbf{u}_1(\mathbf{x}, 0)\psi_1(a)]\psi_0(a)w(a) = 0 \quad \forall(\mathbf{x}, t=0), \tag{12} $$

and

$$ \int_{\Gamma} [\mathbf{u}_0(\mathbf{x}, t)\psi_0(a) + \mathbf{u}_1(\mathbf{x}, t)\psi_1(a)]\psi_0(a)w(a) = 0 \quad \forall(\mathbf{x}, t) \in \Gamma. \tag{13} $$

Hence, Equation (12) gives that for all $\mathbf{x} \in \Omega$ and $t = 0$, $\mathbf{u}_0(\mathbf{x}, 0) = 0$; and Equation (13) gives that for all $\mathbf{x} \in \Gamma$ and $t \ge 0$, $\mathbf{u}_0(\mathbf{x}, t) = 0$. All these equations are now grouped into the following system

$$ \left\{ \begin{array}{ll} \partial_t \mathbf{u}_0(\mathbf{x}, t) - v \Delta \mathbf{u}_0(\mathbf{x}, t) + (\mathbf{u}_0(\mathbf{x}, t) \cdot \nabla) \mathbf{u}_0(\mathbf{x}, t) & \\ + (\mathbf{u}_1(\mathbf{x}, t) \cdot \nabla) \mathbf{u}_1(\mathbf{x}, t) + \frac{1}{\rho} \nabla p_0(\mathbf{x}, t) = \mathbf{f}_0 & \Omega, t > 0, \\ \nabla \cdot \mathbf{u}_0(\mathbf{x}, t) = \mathbf{0} & \Omega, t > 0, \\ \mathbf{u}_0(\mathbf{x}, 0) = \mathbf{0} & \Omega, t = 0, \\ \mathbf{u}_0(\mathbf{x}, t) = \mathbf{0} & \Gamma, t > 0. \end{array} \right. \tag{14} $$

The weak formulation of Equation (14) is written as follows

$$ \left\{ \begin{array}{ll} (\partial_t \mathbf{u}_0, \mathbf{v}_0) + a(\mathbf{u}_0, \mathbf{v}_0) + \frac{1}{\rho} b(\mathbf{v}_0, p_0) + t(\mathbf{u}_0, \mathbf{u}_0, \mathbf{v}_0) + t(\mathbf{u}_1, \mathbf{u}_1, \mathbf{v}_0) = l_0(\mathbf{v}_0) & \forall \mathbf{v}_0 \in \mathbf{X}, \\ b(\mathbf{u}_0, q_0) = 0 & \forall q_0 \in \mathbf{W}. \end{array} \right. \tag{15} $$
For the sake of having a more compact notation, we define the linear term on X as follows

$$l_i(\mathbf{v}) = \int_{\Omega} \mathbf{f}_i \cdot \mathbf{v} \, d\mathbf{x}.$$

The weak formulation (15) can be rewritten in a way which is more convenient for our purposes, in particular for the discretization step. According to the Proposition 1, the trilinear term $t(\mathbf{u}_1, \mathbf{u}_1, \mathbf{v}_0)$ is antisymmetric and equal to $\frac{1}{2}t(\mathbf{u}_1, \mathbf{u}_1, \mathbf{v}_0) - \frac{1}{2}t(\mathbf{u}_1, \mathbf{v}_0, \mathbf{u}_1)$. Hence, by replacing it in the system (15), the following weak formulation is obtained

$$ \left\{ \begin{array}{ll} (\partial_t \mathbf{u}_0, \mathbf{v}_0) + a(\mathbf{u}_0, \mathbf{v}_0) + \frac{1}{\rho} b(\mathbf{v}_0, p_0) + t(\mathbf{u}_0, \mathbf{u}_0, \mathbf{v}_0) + \frac{1}{2} t(\mathbf{u}_1, \mathbf{u}_1, \mathbf{v}_0) & \forall \mathbf{v}_0 \in X, \\ -\frac{1}{2} t(\mathbf{u}_1, \mathbf{v}_0, \mathbf{u}_1) = l_0(\mathbf{v}_0) & \forall q_0 \in W. \\ b(\mathbf{u}_0, q_0) = 0 & \end{array} \right. \quad (16) $$

One can do the same for $i = 1$, and the obtained equations are grouped in the following system

$$ \left\{ \begin{array}{ll} \partial_t \mathbf{u}_1(\mathbf{x}, t) - \nu \Delta \mathbf{u}_1(\mathbf{x}, t) + (\mathbf{u}_0(\mathbf{x}, t) \cdot \nabla) \mathbf{u}_1(\mathbf{x}, t) & \Omega, t > 0, \\ + (\mathbf{u}_1(\mathbf{x}, t) \cdot \nabla) \mathbf{u}_0(\mathbf{x}, t) + \frac{1}{\rho} \nabla p_1(\mathbf{x}, t) = \mathbf{f}_1 & \Omega, t > 0, \\ \nabla \cdot \mathbf{u}_1(\mathbf{x}, t) = \mathbf{0} & \Omega, t > 0, \\ \mathbf{u}_1(\mathbf{x}, 0) = \mathbf{0} & \Omega, t = 0, \\ \mathbf{u}_1(\mathbf{x}, t) = \mathbf{0} & \Gamma, t > 0. \end{array} \right. \quad (17) $$

The weak formulation of Equation (17) is written as follows

$$ \left\{ \begin{array}{ll} (\partial_t \mathbf{u}_1, \mathbf{v}_1) + a(\mathbf{u}_1, \mathbf{v}_1) + \frac{1}{\rho} b(\mathbf{v}_1, p_1) + t(\mathbf{u}_0, \mathbf{u}_1, \mathbf{v}_1) + t(\mathbf{u}_1, \mathbf{u}_0, \mathbf{v}_1) = l_1(\mathbf{v}_1) & \forall \mathbf{v}_1 \in X, \\ b(\mathbf{u}_1, q_1) = 0 & \forall q_1 \in W. \end{array} \right. \quad (18) $$

One again, the trilinear term $t(\mathbf{u}_1, \mathbf{u}_0, \mathbf{v}_1)$ is equal to $\frac{1}{2}t(\mathbf{u}_1, \mathbf{u}_0, \mathbf{v}_1) - \frac{1}{2}t(\mathbf{u}_1, \mathbf{v}_1, \mathbf{u}_0)$, by replacing it in (18) the following weak formulation is obtained

$$ \left\{ \begin{array}{ll} (\partial_t \mathbf{u}_1, \mathbf{v}_1) + a(\mathbf{u}_1, \mathbf{v}_1) + \frac{1}{\rho} b(\mathbf{v}_1, p_1) + t(\mathbf{u}_0, \mathbf{u}_1, \mathbf{v}_1) + \frac{1}{2} t(\mathbf{u}_1, \mathbf{u}_0, \mathbf{v}_1) & \forall \mathbf{v}_1 \in X, \\ -\frac{1}{2} t(\mathbf{u}_1, \mathbf{v}_1, \mathbf{u}_0) = l_1(\mathbf{v}_1) & \forall q_1 \in W. \\ b(\mathbf{u}_1, q_1) = 0 & \end{array} \right. \quad (19) $$

The two formulations (16) and (19) are written because they will be used in Section 6 to discretize the sensitivity of the Navier-Stokes equations according to the FEV scheme and to study the stability of the discrete sensitivity of the Navier-Stokes equations.

One can notice that the systems (14) and (17) are coupled. Therefore, this makes them difficult to solve numerically. To simplify the resolution of these two systems, we perform a decoupling, that is, we find two new systems of PDEs, such that the first system is independent of the second one. The second system will depend on the solution of the first; therefore, we obtain a triangular system. This is done using the following change of variables

$$ \mathbf{u}_1 = \sigma \tilde{\mathbf{u}}_1, \quad p_1 = \sigma \tilde{p}_1, \quad \text{and} \quad \mathbf{f}_1 = \sigma \tilde{\mathbf{f}}_1, $$

then these quantities are inserted in system (17). This new system (17) is subtracted from system (14) and according to the following change of variables

$$ \mathbf{u}_0 = \tilde{\mathbf{u}}_0 + \sigma \tilde{\mathbf{u}}_1, \quad p_0 = \tilde{p}_0 + \sigma \tilde{p}_1, \quad \text{and} \quad \mathbf{f}_0 = \tilde{\mathbf{f}}_0 + \sigma \tilde{\mathbf{f}}_1. $$
Finally, two decoupled systems (20) and (21) equivalent to the systems (14) and (17) are obtained

$$
\left\{
\begin{array}{ll}
\displaystyle \partial_t \tilde{u}_0(\mathbf{x}, t) - \nu \Delta \tilde{u}_0(\mathbf{x}, t) + (\tilde{u}_0(\mathbf{x}, t) \cdot \nabla) \tilde{u}_0(\mathbf{x}, t) + \frac{1}{\rho} \nabla \tilde{p}_0(\mathbf{x}, t) = \tilde{f}_0 & \Omega, t > 0, \\
\displaystyle \nabla \cdot \tilde{u}_0(\mathbf{x}, t) = \mathbf{0} & \Omega, t > 0, \\
\displaystyle \tilde{u}_0(\mathbf{x}, 0) = \mathbf{0} & \Omega, t = 0, \\
\displaystyle \tilde{u}_0(\mathbf{x}, t) = \mathbf{0} & \Gamma, t > 0,
\end{array}
\right.
\tag{20}
$$

with $\tilde{u}_0$, $\tilde{p}_0$ and $\tilde{f}_0$ the velocity, pressure and external force respectively.

$$
\left\{
\begin{array}{ll}
\displaystyle \partial_t \tilde{u}_1(\mathbf{x}, t) - \nu \Delta \tilde{u}_1(\mathbf{x}, t) + (\tilde{u}_0(\mathbf{x}, t) \cdot \nabla) \tilde{u}_1(\mathbf{x}, t) + (\tilde{u}_1(\mathbf{x}, t) \cdot \nabla) \tilde{u}_0(\mathbf{x}, t) & \\
\displaystyle \qquad + 2\sigma(\tilde{u}_1(\mathbf{x}, t) \cdot \nabla) \tilde{u}_1(\mathbf{x}, t) + \frac{1}{\rho} \nabla \tilde{p}_1 = \tilde{f}_1 & \Omega, t > 0, \\
\displaystyle \nabla \cdot \tilde{u}_1(\mathbf{x}, t) = \mathbf{0} & \Omega, t > 0, \\
\displaystyle \tilde{u}_1(\mathbf{x}, 0) = \mathbf{0} & \Omega, t = 0, \\
\displaystyle \tilde{u}_1(\mathbf{x}, t) = \mathbf{0} & \Gamma, t > 0.
\end{array}
\right.
\tag{21}
$$

Remark 2. The system (20) corresponds to the Navier–Stokes equations (1). The system (21) represents the first-order sensitivity of the Navier–Stokes equations, with $\tilde{u}_1$, $\tilde{p}_1$, and $\tilde{f}_1$ the sensitivity of the velocity, pressure, and external force respectively. These equations will also be referred to as sensitivity equations in the following. The system (21) can be solved once system (20) is solved and the value of $\tilde{u}_0$ is computed.

A stability estimate of the solutions $\mathbf{u}_0$ and $\mathbf{u}_1$ is established in the following proposition.

**Proposition 2.** Let $\mathbf{f}_0$ and $\mathbf{f}_1$ be stationary, the following stability estimate holds

$$
\|\mathbf{u}_0\|^2 + \|\mathbf{u}_1\|^2 + \nu \int_0^t \|\nabla \mathbf{u}_0(s)\|^2 + \|\nabla \mathbf{u}_1(s)\|^2 ds \le \min(k_1(t), k_2(t)) (\|\mathbf{f}_0\|^2 + \|\mathbf{f}_1\|^2), \quad (22)
$$

where $\| \cdot \|$ is the $L^2$ norm in space, $k_1(t) = e^t - 1$, $k_2(t) = \frac{4c^2}{\nu}t$, and $c = \frac{\mathrm{Diam}(\Omega)}{\pi}$ the Poincaré constant $^{28,29}$.

*Proof.* The stability estimate, Equation (22), is obtained by replacing $\mathbf{v}_0$ by $\mathbf{u}_0$ and $q_0$ by $p_0$ in Equation (15) and by replacing $\mathbf{v}_1$ by $\mathbf{u}_1$ and $q_1$ by $p_1$ in Equation (18), obtaining

$$
(\partial_t \mathbf{u}_0, \mathbf{u}_0) = \int_{\Omega} \partial_t \frac{\mathbf{u}_0^2}{2} = \frac{1}{2} \frac{d}{dt} \|\mathbf{u}_0\|^2, \quad a(\mathbf{u}_0, \mathbf{u}_0) = -\nu \int_{\Omega} \mathbf{u}_0 \Delta \mathbf{u}_0 = \nu \|\nabla \mathbf{u}_0\|^2, \quad (23)
$$

$$
(\partial_t \mathbf{u}_1, \mathbf{u}_1) = \int_{\Omega} \partial_t \frac{\mathbf{u}_1^2}{2} = \frac{1}{2} \frac{d}{dt} \|\mathbf{u}_1\|^2, \quad a(\mathbf{u}_1, \mathbf{u}_1) = \nu \|\nabla \mathbf{u}_1\|^2, \quad (24)
$$

and thanks to Proposition 1, $t(\mathbf{u}_0, \mathbf{u}_0, \mathbf{u}_0) = t(\mathbf{u}_1, \mathbf{u}_1, \mathbf{u}_1) = 0$.

Using Cauchy–Schwarz and Young inequalities, the linear forms $l_0(\mathbf{u}_0)$ and $l_1(\mathbf{u}_1)$ are bounded as follows

$$
l_0(\mathbf{u}_0) = \int_{\Omega} \mathbf{f}_0 \cdot \mathbf{u}_0 \le \|\mathbf{f}_0\| \|\mathbf{u}_0\| \le \frac{1}{2} \|\mathbf{f}_0\|^2 + \frac{1}{2} \|\mathbf{u}_0\|^2,
$$

$$
l_1(\mathbf{u}_1) = \int_{\Omega} \mathbf{f}_1 \cdot \mathbf{u}_1 \le \|\mathbf{f}_1\| \|\mathbf{u}_1\| \le \frac{1}{2} \|\mathbf{f}_1\|^2 + \frac{1}{2} \|\mathbf{u}_1\|^2.
$$

According to Equation (15) and Equation (18) the following inequalities hold

$$
\frac{1}{2} \frac{d}{dt} \|\mathbf{u}_0\|^2 + \nu \|\nabla \mathbf{u}_0\|^2 + t(\mathbf{u}_1, \mathbf{u}_1, \mathbf{u}_0) \le \frac{1}{2} \|\mathbf{f}_0\|^2 + \frac{1}{2} \|\mathbf{u}_0\|^2, \quad (25)
$$

$$
\frac{1}{2} \frac{d}{dt} \|\mathbf{u}_1\|^2 + \nu \|\nabla \mathbf{u}_1\|^2 + t(\mathbf{u}_1, \mathbf{u}_0, \mathbf{u}_1) \le \frac{1}{2} \|\mathbf{f}_1\|^2 + \frac{1}{2} \|\mathbf{u}_1\|^2. \quad (26)
$$
By summing Equation (25) and Equation (26) and multiplying by 2, knowing that and $t(\mathbf{u}_1, \mathbf{u}_1, \mathbf{u}_0) = -t(\mathbf{u}_1, \mathbf{u}_0, \mathbf{u}_1)$, one has

$$ \begin{aligned} & \frac{d}{dt} (\|\mathbf{u}_0\|^2 + \|\mathbf{u}_1\|^2) + 2\nu(\|\nabla \mathbf{u}_0\|^2 + \|\nabla \mathbf{u}_1\|^2) \\ & \leq \|\mathbf{u}_0\|^2 + \|\mathbf{u}_1\|^2 + \|\mathbf{f}_0\|^2 + \|\mathbf{f}_1\|^2. \end{aligned} \quad (27) $$

By integrating Equation (27) in time, one has

$$ \begin{aligned} & \|\mathbf{u}_0(t)\|^2 + \|\mathbf{u}_1(t)\|^2 + 2\nu \int_0^t \|\nabla \mathbf{u}_0(s)\|^2 + \|\nabla \mathbf{u}_1(s)\|^2 ds \\ & \leq \int_0^t \|\mathbf{u}_0(s)\|^2 + \|\mathbf{u}_1(s)\|^2 ds + \int_0^t \|\mathbf{f}_0\|^2 + \|\mathbf{f}_1\|^2 ds. \end{aligned} \quad (28) $$

Knowing that $\mathbf{f}_0$ and $\mathbf{f}_1$ are stationary and by setting $z(t) = \int_0^t \|\mathbf{u}_0(s)\|^2 + \|\mathbf{u}_1(s)\|^2 ds$ then replacing it in Equation (28) and multiplying by $e^{-t}$, one has

$$ e^{-t}(z'(t) - z(t)) \leq e^{-t}t(\|\mathbf{f}_0\|^2 + \|\mathbf{f}_1\|^2). \quad (29) $$

Integrating Equation (29) in time and multiplying by $e^t$ gives

$$ z(t) \leq (e^t - 1 - t)(\|\mathbf{f}_0\|^2 + \|\mathbf{f}_1\|^2). \quad (30) $$

Finally, $z(t)$ is replaced in Equation (28)

$$ \begin{aligned} & \|\mathbf{u}_0(t)\|^2 + \|\mathbf{u}_1(t)\|^2 + \nu \int_0^t \|\nabla \mathbf{u}_0(s)\|^2 + \|\nabla \mathbf{u}_1(s)\|^2 ds \\ & \leq k_1(t)(\|\mathbf{f}_0\|^2 + \|\mathbf{f}_1\|^2), \end{aligned} \quad (31) $$

with $k_1(t) = e^t - 1$.

According to Cauchy-Schwarz and Poincaré inequalities, the linear forms $l_0(\mathbf{u}_0)$ and $l_1(\mathbf{u}_1)$ can also be bounded as follows

$$ l_0(\mathbf{u}_0) = \int_{\Omega} \mathbf{f}_0 \cdot \mathbf{u}_0 \leq c \|\mathbf{f}_0\| \|\nabla \mathbf{u}_0\| \leq \frac{2c^2}{\nu} \|\mathbf{f}_0\|^2 + \frac{\nu}{2} \|\nabla \mathbf{u}_0\|^2, \quad (32) $$

$$ l_1(\mathbf{u}_1) = \int_{\Omega} \mathbf{f}_1 \cdot \mathbf{u}_1 \leq c \|\mathbf{f}_1\| \|\nabla \mathbf{u}_1\| \leq \frac{2c^2}{\nu} \|\mathbf{f}_1\|^2 + \frac{\nu}{2} \|\nabla \mathbf{u}_1\|^2. \quad (33) $$

By plugging the inequalities (23) and (32) in Equations (15) and (24), (33) in Equation (18) the following inequalities hold

$$ \frac{1}{2} \frac{d}{dt} \|\mathbf{u}_0\|^2 + \nu \|\nabla \mathbf{u}_0\|^2 \leq \frac{2c^2}{\nu} \|\mathbf{f}_0\|^2 + \frac{\nu}{2} \|\nabla \mathbf{u}_0\|^2, \quad (34) $$

$$ \frac{1}{2} \frac{d}{dt} \|\mathbf{u}_1\|^2 + \nu \|\nabla \mathbf{u}_1\|^2 \leq \frac{2c^2}{\nu} \|\mathbf{f}_1\|^2 + \frac{\nu}{2} \|\nabla \mathbf{u}_1\|^2. \quad (35) $$

By summing Equation (34) and Equation (35) one has

$$ \begin{aligned} & \frac{1}{2} \frac{d}{dt} (\|\mathbf{u}_0\|^2 + \|\mathbf{u}_1\|^2) + \nu (\|\nabla \mathbf{u}_0\|^2 + \|\nabla \mathbf{u}_1\|^2) \leq \frac{2c^2}{\nu} \|\mathbf{f}_0\|^2 + \frac{\nu}{2} \|\nabla \mathbf{u}_0\|^2 \\ & \qquad + \frac{2c^2}{\nu} \|\mathbf{f}_1\|^2 + \frac{\nu}{2} \|\nabla \mathbf{u}_1\|^2. \end{aligned} \quad (36) $$
Integrating Equation (36) in time and multiplying by 2 gives

$$ \|\mathbf{u}_0(t)\|^2 + \|\mathbf{u}_1(t)\|^2 + \nu \int_0^t \|\nabla \mathbf{u}_0(s)\|^2 + \|\nabla \mathbf{u}_1(s)\|^2 ds \le k_2(t)(\|\mathbf{f}_0\|^2 + \|\mathbf{f}_1\|^2), \quad (37) $$

with $k_2(t) = t^{2c^2/\nu}$. Finally, according to the Equation (31) and Equation (37) one has

$$ \|\mathbf{u}_0\|^2 + \|\mathbf{u}_1\|^2 + \nu \int_0^t \|\nabla \mathbf{u}_0(s)\|^2 + \|\nabla \mathbf{u}_1(s)\|^2 ds \le \min(k_1(t), k_2(t))(\|\mathbf{f}_0\|^2 + \|\mathbf{f}_1\|^2). \quad \blacksquare $$

# 4 | DESCRIPTION OF THE FINITE ELEMENT-VOLUME METHOD

In this section, we describe the FEV method that we will use to discretize the sensitivity of the Navier-Stokes equations. This discretization is implemented in the industrial code TrioCFD.20,21,30

TrioCFD is an open-source software developed by the CEA (French Alternative Energies and Atomic Energy Commission). It is designed to tackle a wide range of complex challenges, such as turbulent fluid dynamics, fluid-solid coupling, and multiphase flows and so forth. This software addresses problems in the nuclear energy sector through mathematical and computational modeling. FEV is one of the most commonly used spatial discretizations in TrioCFD. The FEV was first proposed by Emonot in his thesis,¹⁶ then extended by in References 17 and 18. The FEV method has been extensively applied to Stokes and Navier-Stokes equations and, to a lesser extent, to sensitivity equations. Notably, this scheme needs more numerical analysis concerning the sensitivity of the Navier-Stokes equations. In this section, the FEV method is described in detail. First, the mesh and localization of the degrees of freedom are presented. Then, the different control volumes are described. Finally, the decomposition of velocity and pressure is explained.

## 4.1 | Mesh

Let $h$ be the mesh size and $T_h$ be a regular triangulation of $\Omega$, whose elements are noted $K_l$, in the usual sense where:

* for any $h$, the intersection of two different elements of $T_h$, if not empty, is either a vertex or an edge common to these two elements.

* the ratio between the diameter $h_{K_l}$ of an element $K_l \in T_h$ and the diameter of the inscribed circle is bounded by a constant $r$ independent of $h$.

$N_K$, $N_S$, $N_{FY}$, and $N_E$ denote respectively, the number of elements, vertices, faces of $T_h$, and the elements that have $s_j$ as a vertex. The set of faces is denoted by $F_h$. The vertices are $s_j$, the middle points of faces by $x_i$, the faces by $f_i$, and the centers of mass by $c_l$.

## 4.2 | Localization of the degrees of freedom

Let $P^0(K_l)$ be the space of constant functions on the triangle $K_l$ and $P^1(K_l)$ be the space of linear functions on $K_l$. The two finite element spaces $X_h$ and $W_h$ are defined as the approximation spaces for the velocity $u_h$ and the pressure $p_h$, respectively.

These spaces are introduced as follows

$$ \begin{align*} X_h &= \{v_h \text{ continuous at } \mathbf{x}_i : \forall K_l \in T_h, v_h|_{K_l} \in P_1(K_l)\}, \\ \mathbf{X}_h &= \{\mathbf{v}_h = (v_x, v_y) : v_x, v_y \in V_h\}, \\ W_h &= \{p_h : \forall K_l \in T_h, p_h|_{K_l} \in P^0(K_l) \oplus P^1(K_l)\}, \\ \text{and, } W_{h,0} &= W_h \cap L_0^2(\Omega), \text{ with } L_0^2(\Omega) = \{q_h \in L^2(\Omega); \int_{\Omega} q_h \, d\mathbf{x}\}. \end{align*} $$
The space $\mathbf{X}_h$ is endowed with the following semi-norm $|\mathbf{u}_h|_{h,\Omega} = (\sum_{l=1}^{N_K} |\mathbf{u}_h|^2_{H^1(K_l)})^{\frac{1}{2}}$.

The velocity degrees of freedom are located on the middle points of the faces of the triangles. The pressure degrees of freedom are located at the vertices and the center of gravity of the elements, see Figure 1. The space $\mathbf{X}_h$ is therefore not included in $\mathbf{X} = H_0^1(\Omega)^2$, making this non-conforming.

Let $\phi_i$ be the one-dimensional basis function of $\mathbf{X}_h$ which verifies $\phi_i(\mathbf{x}_j) = \delta_{ij}$. On each element $K_l$, the basis functions are written, $\phi_i|_{K_l} = 1 - 2\lambda_i|_{K_l}$ where $\lambda_i|_{K_l}$ is the barycentric coordinate of $K_l$ associated to the vertex facing the node $\mathbf{x}_i$.

## 4.3 | Description of the control volumes

To apply the FEV method, a mesh of the domain and a finite element space are needed. The solution is characterized by its values at the nodes, and a control volume is associated with each node.

For the Navier-Stokes equations, the control volume associated to each velocity node located at the center of a face is constructed by joining the centers of gravity of the elements having a face in common with the vertices of that face, see Figure 2B. The $\kappa_i$ in each element is defined by the following relation (see Figure 2A):

$$ \kappa_i = \{ \mathbf{x} \in K : \phi_i(\mathbf{x}) \ge \phi_j(\mathbf{x}) \forall i \ne j \}. $$

A control volume $\omega_i$ is defined as $\omega_i = \bigcup_{j_i \in K_i} \kappa_i$.

The two types of control volumes for the mass conservation equation are

* the first one is associated with pressure node at the gravity center $c_l$ of an element $K_l$ is $K_l$.

* the second one, denoted by $\Pi_{s_j}$, is associated with the pressure node in $s_j$.

Different choices can be made for the construction of control volumes. These choices and their consequences on the definition of the bilinear forms are detailed in Reference 18. The control volume $\Pi_{s_j}$ is defined in the following way:

<figure><img src="image_3.png" /><figcaption>FIGURE 1 Localization of the degrees of freedom for the element $P_{NC}^1 \setminus P^0 + P^1$. [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 2 The figure (A) represents the $\kappa_i$ attached to the node $\mathbf{x}_i$ and the figure (B) represents the control volume $\omega_i$ (the dotted area) also attached to $\mathbf{x}_i$. [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIGURE 3 The control volume K associated with the pressure node on c_l (the hatched area) and the control volume Π_s_j associated with the vertex (the dotted area). [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

$\Pi_{s_j} = \bigcup_{K_l \in T_h} \{x \in K_l; \lambda_{s_j}(x) \ge \frac{1}{2}\}$. This control volume $\Pi_{s_j}$ is constructed by joining on all the elements containing $s_j$ the centers of the faces of which $s_j$ is a vertex, see Figure 3.

## 4.4 Decomposition of velocity and pressure spaces

The two spaces $X_h$ and $W_h$ are each decomposed into two subspaces as follows:

* $X_h = X_h^T + X_h^N$ with $X_h^T$ the space of velocities $P_{NC}^1$ tangent to the faces and $X_h^N$ the space of velocities $P_{NC}^1$ normal to the faces.

* $W_h = W_h^0 + W_h^1$ with $W_h^0$ the space of pressures spanned by $P^0(K_l)$ and $W_h^1$ the space of pressures spanned by $P^1(K_l)$.

Any function $u_h$ of $X_h$ can be written

$$ \forall \mathbf{x} \in \Omega, \quad \mathbf{u}_h(\mathbf{x}) = \begin{bmatrix} \sum_{i=1}^{N_F} u_h^x(\mathbf{x}_i)\phi_i(\mathbf{x}) \\ \sum_{i=1}^{N_F} u_h^y(\mathbf{x}_i)\phi_i(\mathbf{x}) \end{bmatrix}. \qquad (38) $$

Any function $p_h$ of $W_h$ is written

$$ \forall \mathbf{x} \in \Omega, \quad p_h(\mathbf{x}) = \sum_{j=1}^{N_S} p_h(\mathbf{s}_j)\lambda_{\mathbf{s}_j}(\mathbf{x}) + \sum_{l=1}^{N_K} \left( p_h(\mathbf{c}_l) - \frac{1}{3} \sum_{i=1}^{N_S} p_h(\mathbf{s}_j) \right) \mathbf{1}_{K_l}(\mathbf{x}). \qquad (39) $$

As $W_h = W_h^0 + W_h^1$, $p_h$ is written as follows

$$ \forall \mathbf{x} \in \Omega, \quad p_h(\mathbf{x}) = p_h^1(\mathbf{x}) + p_h^0(\mathbf{x}), $$

with $p_h^1 \in W_h^1$ and $p_h^0 \in W_h^0$ defined in the following way

$$ p_h^1 = \sum_{i=1}^{N_S} p_h(\mathbf{s}_j)\lambda_{\mathbf{s}_j} = \sum_{i=1}^{N_S} p_h^i \lambda_{\mathbf{s}_j}, $$

$$ p_h^0 = \sum_{l=1}^{N_K} \left( p_h(\mathbf{c}_l) - \frac{1}{3} \sum_{j=1}^{N_S} p_h(\mathbf{s}_j) \right) \mathbf{1}_{K_l} = \sum_{l=1}^{N_K} p_h^{K_l} \mathbf{1}_{K_l}. $$
The notations introduced are

$$ \forall j = 1, \dots, N_s, \quad p_h^{\mathbf{s}_j} = p_h(\mathbf{s}_j), $$

$$ \forall l = 1, \dots, N_{K_l}, \quad p_h^{K_l} = p_h(\mathbf{c}_l) - \frac{1}{3} \sum_{j=1}^{N_s} p_h(\mathbf{s}_j). $$

The FEV method is a widely employed discretization technique for PDEs. A contributing factor to its popularity is its capacity to maintain local mass conservation, a vital aspect in tackling particular physical problems. The method can be regarded as a modification of the finite element method, except for the source term,¹⁶ when used to solve the Stokes equations, and it makes use of established theories from the finite element method, especially to show stability and the convergence of the FEV scheme for the Stokes equations.¹⁸

# 5 | DISCRETE NAVIER-STOKES EQUATIONS

In this section, the Navier-Stokes equations, Equation (1), are discretized using the FEV method, and a stability estimate is established for the discrete Navier-Stokes equations. These results will be generalized to the sensitivity equations.

## 5.1 | Spatial discretization

In order to have a discrete formulation, the conservation of momentum equation is integrated over $\omega_i$ and then multiplied by a test function associated with the node $\mathbf{x}_i$. The conservation of mass equation is integrated over $K_l$ then multiplied by a test function associated with $\mathbf{c}_l$; it is also integrated over $\Pi_s$, then multiplied by a test function associated with the node $\mathbf{s}_j$. Finally, we sum over all nodes. The discrete Navier-Stokes problem formulated in FEV is written as follows: find $(u_h, p_h) \in \mathbf{X}_h \times W_{h,0}$ such that

$$ \left\{ \begin{array}{ll} a_h^V(\mathbf{u}_h, \mathbf{v}_h) + \frac{1}{\rho} b_h^V(\mathbf{v}_h, p_h) + t_h^V(\mathbf{u}_h, \mathbf{u}_h, \mathbf{v}_h) = l_h^V(\mathbf{v}_h) & \forall \mathbf{v}_h \in \mathbf{X}_h, \\ c_h^V(\mathbf{u}_h, q_h) = 0 & \forall q_h \in W_{h,0}. \end{array} \right. \quad (40) $$

with $a_h^V(\mathbf{u}_h, \mathbf{v}_h)$ the following bilinear form on $\mathbf{X}_h \times \mathbf{X}_h$

$$ a_h^V(\mathbf{u}_h, \mathbf{v}_h) = -\nu \sum_{i=1}^{N_F} \left( \mathbf{v}_h(\mathbf{x}_i) \cdot \int_{\partial\omega_i} \nabla \mathbf{u}_h \cdot \mathbf{n} \, ds \right), $$

$b_h^V(\mathbf{v}_h, p_h)$ and $c_h^V(\mathbf{u}_h, q_h)$ the two bilinear forms on $\mathbf{X}_h \times W_h$

$$ b_h^V(\mathbf{v}_h, p_h) = \sum_{i=1}^{N_F} \left( \mathbf{v}_h(\mathbf{x}_i) \cdot \int_{\partial\omega_i} p_h \mathbf{n} \, ds \right), $$

$$ c_h^V(\mathbf{u}_h, q_h) = \sum_{j=1}^{N_s} q_h^{\mathbf{s}_j} \left( \int_{\partial\Pi_{s_j}} \mathbf{u}_h \cdot \mathbf{n} \, ds + \alpha_1 \sum_{l=1}^{N_F} \int_{\partial K_l} \mathbf{u}_h \cdot \mathbf{n} \, ds \right) + \alpha_2 \sum_{l=1}^{N_K} q_h^{K_l} \left( \int_{\partial K_l} \mathbf{u}_h \cdot \mathbf{n} \, ds \right), $$

$l_h^V(\mathbf{v}_h)$ the linear form on $\mathbf{X}_h$

$$ l_h^V(\mathbf{v}_h) = \sum_{i=1}^{N_F} \left( \mathbf{v}_h(\mathbf{x}_i) \cdot \int_{\omega_i} \mathbf{f} \, dx \right), $$

and $t_h^V(\mathbf{u}_h, \mathbf{v}_h, \mathbf{w}_h)$ the trilinear form on $\mathbf{X}_h \times \mathbf{X}_h \times \mathbf{X}_h$

$$ t_h^V(\mathbf{u}_h, \mathbf{v}_h, \mathbf{w}_h) = \sum_{i=1}^{N_F} \left( \mathbf{w}_h(\mathbf{x}_i) \cdot \int_{\omega_i} (\mathbf{u}_h(\mathbf{x}) \cdot \nabla) \mathbf{v}_h(\mathbf{x}) \, dx \right) = \sum_{i=1}^{N_F} \left( \mathbf{w}_h(\mathbf{x}_i) \cdot \int_{\partial\omega_i} \mathbf{v}_h(\mathbf{x}) (\mathbf{u}_h(\mathbf{x}) \cdot \mathbf{n}_i(\mathbf{x})) \, ds \right). $$
The two coefficients $\alpha_1$ and $\alpha_2$ are used to make the form $c_h^V(\cdot, \cdot)$ proportional to $b_h^V(\cdot, \cdot)$. For our choice of control volume, the coefficients are $\alpha_1 = 1/4$ and $\alpha_2 = 3/2$. These two coefficients depend on the space dimension and the choice of the mass control volumes; see Reference 18 for the computation of the different options for $\alpha_1$ and $\alpha_2$.

The linear and bilinear forms of the weak formulation (40) are compared in Reference 18 with those of the following discretization of the Navier-Stokes equations, obtained with the finite element method: find $(\mathbf{u}_h, p_h) \in \mathbf{X}_h \times W_{h,0}$ such that

$$ \left\{ \begin{array}{ll} a_h(\mathbf{u}_h, \mathbf{v}_h) + \frac{1}{\rho} b_h(\mathbf{v}_h, p_h) + t_h(\mathbf{u}_h, \mathbf{u}_h, \mathbf{v}_h) = l_h(\mathbf{v}_h) & \forall \mathbf{v}_h \in \mathbf{X}_h, \\ b_h(\mathbf{u}_h, q_h) = 0 & \forall q_h \in W_{h,0}, \end{array} \right. \quad (41) $$

where the bilinear form $a_h(\cdot, \cdot)$, the linear form $l_h(\cdot)$, and the trilinear form $t_h(\cdot, \cdot, \cdot)$ are defined as follows

$$ \begin{aligned} a_h(\mathbf{u}_h, \mathbf{v}_h) &= v \sum_{l=1}^{N_K} \int_{K_l} \nabla \mathbf{u}_h : \nabla \mathbf{v}_h \,dx, \quad l_h(\mathbf{v}_h) = \sum_{l=1}^{N_K} \int_{K_l} \mathbf{f}_h \cdot \mathbf{v}_h \,dx, \\ t_h(\mathbf{u}_h, \mathbf{v}_h, \mathbf{w}_h) &= \sum_{l=1}^{N_K} \int_{K_l} [\mathbf{u}_h \cdot \nabla \mathbf{v}_h] \cdot \mathbf{w}_h \,dx. \end{aligned} $$

To define the bilinear form $b_h(\cdot, \cdot)$ another decomposition of the space $W_h$ is used, $W_h = \widetilde{W}_h \oplus W_h^0$ with

$$ \widetilde{W}_h = \{ q_h \in W_h ; \forall K_l \in \mathcal{T}_h, q_{h|_{K_l}} \in L_0^2(K_l) \}. $$

The orthogonal projection operator $\Pi_h^0 : L_0^2(K_l) \to W_h^0$ is introduced:

$$ \forall p \in L^2(\Omega), \forall q_h \in W_h^0 \quad (p - \Pi_h^0 p; q_h) = 0. $$

Then the second orthogonal projection operator $\tilde{\Pi}_h : L_0^2(K_l) \to \widetilde{W}_h$ is defined:

$$ \forall p \in L^2(\Omega), \forall q_h \in \widetilde{W}_h \quad (p - \tilde{\Pi}_h p; q_h) = 0. $$

Finally, the linear form $b_h$ is defined

$$ b_h(\mathbf{v}, p) = \sum_{l=1}^{N_K} \int_{K_l} \mathbf{v} \cdot \nabla p_h^1 \,dx - \sum_{l=1}^{N_K} \int_{K_l} \nabla \cdot \mathbf{v}(\bar{p}_h - \bar{p}_h^1) \,dx. $$

with $p \in L_0^2(\Omega)$, $\bar{p}_h = \Pi_h^0 p$, and $\tilde{p}_h = \tilde{\Pi}_h p = p_h^1 - \bar{p}_h^1$. The bilinear forms $a_h(\cdot, \cdot)$ and $a_h^V(\cdot, \cdot)$ are identical on the space $\mathbf{X}_h$ 16 and $\forall \mathbf{v}_h \in \mathbf{X}_h$ and $p_h \in W_h$, $b_h(\mathbf{v}_h, p_h) = b_h^V(\mathbf{v}_h, p_h)$. 18

We now define the trilinear term $t_{hh}^V(\mathbf{u}_h, \mathbf{v}_h, \mathbf{w}_h)$

$$ t_h^V(\mathbf{u}_h, \mathbf{v}_h, \mathbf{w}_h) \approx t_{hh}^V(\mathbf{u}_h, \mathbf{v}_h, \mathbf{w}_h) \equiv \sum_{l=1}^{N_F} \left( \mathbf{w}_h(\mathbf{x}_l) \cdot \int_{\partial\omega_l} \mathbf{v}_h(\mathbf{x}_*) (\mathbf{u}_h(\mathbf{x}) \cdot \mathbf{n}_l(\mathbf{x})) ds \right), $$

where $\mathbf{x}_*$ is defined as follows

$$ \mathbf{x}_* = \begin{cases} \mathbf{x}_i & \text{if } \mathbf{u}_h \cdot \mathbf{n}_i > 0 \text{ on } \partial\omega_i \cap \partial\omega_j, \\ \mathbf{x}_j & \text{otherwise.} \end{cases} $$

This corresponds to an upwind-type approach. For each edge of $\omega_i$, $\mathbf{u}_h \cdot \mathbf{n}_i$ is computed and, according to its sign, it is multiplied by either $\mathbf{u}_h(\mathbf{x}_i)$ or $\mathbf{u}_h(\mathbf{x}_j)$, where $\omega_j$ is the control volume adjacent to the considered edge. The upwind-type scheme leads to a CFL condition 31 on the time step $\Delta t$ when considering the time discretization. The trilinear term
$t_{hh}^V(\cdot, \cdot, \cdot)$ possesses an advantageous and indispensable property, which will be presented and verified in the following proposition.

**Proposition 3.** Let $\mathbf{u}_h \in \mathbf{X}_h$ such that $\nabla \cdot \mathbf{u}_h = 0$, then

$$\forall \mathbf{v}_h \in \mathbf{X}_h, \quad t_{hh}^V(\mathbf{u}_h, \mathbf{v}_h, \mathbf{v}_h) \geq 0.$$

*Proof.*

$$t_{hh}^V(\mathbf{u}_h, \mathbf{v}_h, \mathbf{v}_h) = \sum_{j \neq i} \sum_i \mathbf{v}_h(\mathbf{x}_i) \cdot \left( \int_{\partial \omega_i \cap \partial \omega_j} \mathbf{v}_h(\mathbf{x}_*) (\mathbf{u}_h(\mathbf{x}) \cdot \mathbf{n}_i(\mathbf{x})) ds \right)$$

Knowing that $\nabla \cdot \mathbf{u}_h = 0$, the following equality is obtained $\forall \beta \in \mathbb{R}$

$$\begin{aligned} t_{hh}^V(\mathbf{u}_h, \mathbf{v}_h, \mathbf{v}_h) &= \sum_{j \neq i} \sum_i \mathbf{v}_h(\mathbf{x}_i) \cdot \left( \int_{\partial \omega_i \cap \partial \omega_j} \mathbf{v}_h(\mathbf{x}_*) (\mathbf{u}_h(\mathbf{x}) \cdot \mathbf{n}_i(\mathbf{x})) ds \right. \\ & \qquad \left. + \beta \mathbf{v}_h(\mathbf{x}_i) \cdot \int_{\omega_i} \nabla \cdot \mathbf{u}_h dx \right) \\ &= \sum_{j \neq i} \sum_i \mathbf{v}_h(\mathbf{x}_i) \cdot \left( \int_{\partial \omega_i \cap \partial \omega_j} \mathbf{v}_h(\mathbf{x}_*) (\mathbf{u}_h(\mathbf{x}) \cdot \mathbf{n}_i(\mathbf{x})) ds \right. \\ & \qquad \left. + \beta \mathbf{v}_h(\mathbf{x}_i) \cdot \int_{\partial \omega_i \cap \partial \omega_j} \mathbf{u}_h(\mathbf{x}) \cdot \mathbf{n}_i(\mathbf{x}) ds \right) \end{aligned}$$

The proof will be done according to 2 cases.

*Case 1: if $\mathbf{u}_h \cdot \mathbf{n}_i > 0$,*

$$\begin{aligned} t_{hh}^V(\mathbf{u}_h, \mathbf{v}_h, \mathbf{v}_h) &= \sum_{j \neq i} \sum_i \mathbf{v}_h(\mathbf{x}_i) \cdot \left( \int_{\partial \omega_i \cap \partial \omega_j} \mathbf{v}_h(\mathbf{x}_i) (\mathbf{u}_h(\mathbf{x}) \cdot \mathbf{n}_i(\mathbf{x})) ds + \beta \mathbf{v}_h(\mathbf{x}_i) \cdot \int_{\omega_i} \nabla \cdot \mathbf{u}_h dx \right) \\ &= (\mathbf{v}_h^2(\mathbf{x}_i) + \beta \mathbf{v}_h^2(\mathbf{x}_i)) \cdot \left( \int_{\partial \omega_i \cap \partial \omega_j} \mathbf{u}_h(\mathbf{x}) \cdot \mathbf{n}_i(\mathbf{x}) ds \right) \\ & \quad + (\mathbf{v}_h(\mathbf{x}_i) \cdot \mathbf{v}_h(\mathbf{x}_j) + \beta \mathbf{v}_h^2(\mathbf{x}_j)) \cdot \left( \int_{\partial \omega_i \cap \partial \omega_j} \mathbf{u}_h(\mathbf{x}) \cdot \mathbf{n}_j(\mathbf{x}) ds \right) \\ &= (\mathbf{v}_h^2(\mathbf{x}_i) + \beta \mathbf{v}_h^2(\mathbf{x}_i) - \mathbf{v}_h(\mathbf{x}_i) \cdot \mathbf{v}_h(\mathbf{x}_j) - \beta \mathbf{v}_h^2(\mathbf{x}_j)) \cdot \left( \int_{\partial \omega_i \cap \partial \omega_j} \mathbf{u}_h(\mathbf{x}) \cdot \mathbf{n}_i(\mathbf{x}) ds \right) \end{aligned}$$

$$\text{Let } \beta = -\frac{1}{2}, \quad t_{hh}^V(\mathbf{u}_h, \mathbf{v}_h, \mathbf{v}_h) = (\mathbf{v}_h^2(\mathbf{x}_i) - \frac{1}{2}\mathbf{v}_h^2(\mathbf{x}_i) - \mathbf{v}_h(\mathbf{x}_i) \cdot \mathbf{v}_h(\mathbf{x}_j) \\ + \frac{1}{2}\mathbf{v}_h^2(\mathbf{x}_j)) \left( \int_{\partial \omega_i \cap \partial \omega_j} \mathbf{u}_h(\mathbf{x}) \cdot \mathbf{n}_i(\mathbf{x}) ds \right) \\ = \frac{1}{2} |\mathbf{v}_h(\mathbf{x}_i) - \mathbf{v}_h(\mathbf{x}_j)|^2 \left( \int_{\partial \omega_i \cap \partial \omega_j} \mathbf{u}_h(\mathbf{x}) \cdot \mathbf{n}_i(\mathbf{x}) ds \right).$$

Therefore $t_{hh}^V(\mathbf{u}_h, \mathbf{v}_h, \mathbf{v}_h) \geq 0$.

*Case 2: if $\mathbf{u}_h \cdot \mathbf{n}_i < 0$,*

$$\begin{aligned} t_{hh}^V(\mathbf{u}_h, \mathbf{v}_h, \mathbf{v}_h) &= \sum_{j \neq i} \sum_i \mathbf{v}_h(\mathbf{x}_i) \cdot \left( \int_{\partial \omega_i \cap \partial \omega_j} \mathbf{v}_h(\mathbf{x}_j) (\mathbf{u}_h(\mathbf{x}) \cdot \mathbf{n}_i(\mathbf{x})) ds + \beta \mathbf{v}_h(\mathbf{x}_i) \cdot \int_{\omega_i} \nabla \cdot \mathbf{u}_h dx \right) \\ &= (\mathbf{v}_h(\mathbf{x}_i) \cdot \mathbf{v}_h(\mathbf{x}_j) + \beta \mathbf{v}_h^2(\mathbf{x}_i)) \left( \int_{\partial \omega_i \cap \partial \omega_j} \mathbf{u}_h(\mathbf{x}) \cdot \mathbf{n}_i(\mathbf{x}) ds \right) \end{aligned}$$
$$ + (\mathbf{v}_h^2(\mathbf{x}_j) + \beta \mathbf{v}_h^2(\mathbf{x}_j)) \left( \int_{\partial \omega_i \cap \partial \omega_j} \mathbf{u}_h(\mathbf{x}) \cdot \mathbf{n}_j(\mathbf{x}) \, ds \right) \\ = (\beta \mathbf{v}_h^2(\mathbf{x}_i) - \beta \mathbf{v}_h^2(\mathbf{x}_j) + \mathbf{v}_h(\mathbf{x}_i)\mathbf{v}_h(\mathbf{x}_j) - \mathbf{v}_h^2(\mathbf{x}_j)) \left( \int_{\partial \omega_i \cap \partial \omega_j} \mathbf{u}_h(\mathbf{x}) \cdot \mathbf{n}_i(\mathbf{x}) \, ds \right). $$

$$ \text{Let } \beta = -\frac{1}{2}, \quad t_{hh}^V(\mathbf{u}_h, \mathbf{v}_h, \mathbf{v}_h) = \left( -\frac{1}{2}\mathbf{v}_h^2(\mathbf{x}_i) + \mathbf{v}_h(\mathbf{x}_i) \cdot \mathbf{v}_h(\mathbf{x}_j) + \frac{1}{2}\mathbf{v}_h^2(\mathbf{x}_j) \right. \\ \left. - \mathbf{v}_h^2(\mathbf{x}_j) \right) \left( \int_{\partial \omega_i \cap \partial \omega_j} \mathbf{u}_h(\mathbf{x}) \cdot \mathbf{n}_i(\mathbf{x}) \, ds \right) \\ = -\frac{1}{2} |\mathbf{v}_h(\mathbf{x}_i) - \mathbf{v}_h(\mathbf{x}_j)|^2 \left( \int_{\partial \omega_i \cap \partial \omega_j} \mathbf{u}_h(\mathbf{x}) \cdot \mathbf{n}_i(\mathbf{x}) \, ds \right). $$

Therefore $t_{hh}^V(\mathbf{u}_h, \mathbf{v}_h, \mathbf{v}_h) \ge 0$.

<figure><img src="image_7.png" /></figure>

## 5.2 | Time discretization and stability estimate

For the time discretization, an Euler scheme is applied.³² Let $t^n$ be the $n^{th}$ time step and $\mathbf{u}_h^n$ be an approximation of $\mathbf{u}(\mathbf{x}, t^n)$. The interval $[0, T]$ is divided into $N$ intervals having a length $\Delta t^n = t^{n+1} - t^n$. By coupling this with the spatial scheme described in the previous subsection, the following system is obtained find $\mathbf{u}_h^{n+1}$ and $p_h^{n+1}$ (when $\mathbf{u}_h^0, \dots, \mathbf{u}_h^n$ and $p_h^0, \dots, p_h^n$ are known) such that:

$$ \left\{ \begin{array}{l} \displaystyle \sum_{i=1}^{N_F} \int_{\omega_i} \frac{\mathbf{u}_h^{n+1} - \mathbf{u}_h^n}{\Delta t^n} \cdot \mathbf{v}_h + a_h^V(\mathbf{u}_h^{n+1}, \mathbf{v}_h) + t_{hh}^V(\mathbf{u}_h^n, \mathbf{u}_h^{n+1}, \mathbf{v}_h) \\ \displaystyle \qquad + \frac{1}{\rho} b_h^V(\mathbf{v}_h, p_h^{n+1}) = l_h^V(\mathbf{v}_h) \qquad \forall \mathbf{v}_h \in \mathbf{X}_h, \\ \\ c_h^V(\mathbf{u}_h^{n+1}, q_h) = 0 \qquad \forall q_h \in W_{h,0}. \end{array} \right. \quad (42) $$

The following proposition provides a stability estimate for the solution $\mathbf{u}_h^{n+1}$. It is achieved by comparing the bilinear forms of the FEV discretized equations with those of the equations discretized according to the finite element method. The bilinear terms $a_h^V(\cdot, \cdot)$ and $b_h^V(\cdot, \cdot)$ are replaced in system (42) respectively with $a_h(\cdot, \cdot)$ and $b_h(\cdot, \cdot)$:

$$ \left\{ \begin{array}{l} \displaystyle \sum_{i=1}^{N_F} \int_{\omega_i} \frac{\mathbf{u}_h^{n+1} - \mathbf{u}_h^n}{\Delta t^n} \cdot \mathbf{v}_h + a_h(\mathbf{u}_h^{n+1}, \mathbf{v}_h) + \frac{1}{\rho} b_h(\mathbf{v}_h, p_h^{n+1}) \\ \displaystyle \qquad + t_{hh}^V(\mathbf{u}_h^n, \mathbf{u}_h^{n+1}, \mathbf{v}_h) = l_h(\mathbf{v}_h) \qquad \forall \mathbf{v}_h \in \mathbf{X}_h, \\ \\ b_h(\mathbf{u}_h^{n+1}, q_h) = 0 \qquad \forall q_h \in W_{h,0}. \end{array} \right. \quad (43) $$

The proof of the following proposition is presented in Appendix A.

**Proposition 4.** Let $\mathbf{f}_h$ be stationary, the following stability estimate holds

$$ \begin{align} \|\mathbf{u}_h^{n+1}\|^2 + \nu \int_0^{t^{n+1}} \|\nabla \mathbf{u}_h(s)\|^2 ds + 2 \int_0^{t^n} t_{hh}^V(\mathbf{u}_h(s), \mathbf{u}_h(s + \Delta t), \mathbf{u}_h(s + \Delta t)) ds \tag{44} \\ \leq \min(k_1(t^{n+1}), k_2(t^{n+1})) \|\mathbf{f}_h\|^2. \nonumber \end{align} $$

with $k_1(t^{n+1}) = 12(e^{t^{n+1}} - 1)$, $k_2(t^{n+1}) = \frac{48c^2}{\nu} t^{n+1}$, and c is the Poincaré constant ²⁸,²⁹.

# 6 | DISCRETE SENSITIVITY OF THE NAVIER-STOKES EQUATIONS

This section considers the sensitivity of the Navier-Stokes equations from Section 3.2. These equations are discretized using the FEV scheme, and a stability estimate of the discrete sensitivity is established.
## 6.1 | Spatial discretization

There are several ways to discretize the sensitivity equations.

**First method:** Equation (14) and Equation (17) are considered and only the spatial terms are treated. They are discretized using the method described in Section 4. The discrete form of the system (14) formulated in FEV is: find $(\mathbf{u}_{0,h}, p_{0,h}) \in \mathbf{X}_h \times W_{h,0}$ such that

$$ egin{cases} a_h^V(\mathbf{u}_{0,h}, \mathbf{v}_{0,h}) + \frac{1}{\rho} b_h^V(\mathbf{v}_{0,h}, p_{0,h}) + t_{hh}^V(\mathbf{u}_{0,h}, \mathbf{u}_{1,h}, \mathbf{v}_{0,h}) \\ + t_{hh}^V(\mathbf{u}_{1,h}, \mathbf{u}_{1,h}, \mathbf{v}_{0,h}) = l_{0,h}^V(\mathbf{v}_{0,h}) & \forall \mathbf{v}_{0,h} \in \mathbf{X}_h, \\ c_h^V(\mathbf{u}_{0,h}, q_{0,h}) = 0 & \forall q_{0,h} \in W_{h,0}. \end{cases} \quad (45) $$

The discretization of the system (17) formulated in FEV is then written: find $(\mathbf{u}_{1,h}, p_{1,h}) \in \mathbf{X}_h \times W_{h,0}$ such that

$$ egin{cases} a_h^V(\mathbf{u}_{1,h}, \mathbf{v}_{1,h}) + \frac{1}{\rho} b_h^V(\mathbf{v}_{1,h}, p_{1,h}) + t_{hh}^V(\mathbf{u}_{0,h}, \mathbf{u}_{1,h}, \mathbf{v}_{1,h}) \\ + t_{hh}^V(\mathbf{u}_{1,h}, \mathbf{u}_{0,h}, \mathbf{v}_{1,h}) = l_{1,h}^V(\mathbf{v}_{1,h}) & \forall \mathbf{v}_{1,h} \in \mathbf{X}_h, \\ c_h^V(\mathbf{u}_{1,h}, q_{1,h}) = 0 & \forall q_{1,h} \in W_{h,0}. \end{cases} \quad (46) $$

The linear term $l_{m,h}^V(\cdot)$ is defined on $\mathbf{X}_h$ as follows

$$ l_{m,h}^V(\mathbf{v}_h) = \sum_{i=1}^{N_F} \int_{\omega_i} \mathbf{f}_{m,h} \cdot \mathbf{v}_h \, d\mathbf{x}, \quad \text{with } m \in \{0,1\}. $$

In Appendix B, we show the continuity of these linear terms.

**Second method:** The two weak formulations (16) and (19) are considered. Only the spatial part of these systems is treated. To discretize them according to the FEV scheme, the linear term $l_m(\cdot)$, the bilinear terms $a(\cdot, \cdot)$ and $b(\cdot, \cdot)$ and the trilinear one $t(\cdot, \cdot, \cdot)$ are respectively replaced by the discrete linear term $l_{m,h}^V(\cdot)$, the bilinear terms $a_h^V(\cdot, \cdot)$ and $b_h^V(\cdot, \cdot)$ and the trilinear one $t_{hh}^V(\cdot, \cdot, \cdot)$. The discrete version of the system (16) is then written: find $(\mathbf{u}_{0,h}, p_{0,h}) \in \mathbf{X}_h \times W_{h,0}$ such that

$$ egin{cases} a_h^V(\mathbf{u}_{0,h}, \mathbf{v}_{0,h}) + \frac{1}{\rho} b_h^V(\mathbf{v}_{0,h}, p_{0,h}) + t_{hh}^V(\mathbf{u}_{0,h}, \mathbf{u}_{1,h}, \mathbf{v}_{0,h}) \\ + \frac{1}{2} t_{hh}^V(\mathbf{u}_{1,h}, \mathbf{u}_{1,h}, \mathbf{v}_{0,h}) - \frac{1}{2} t_{hh}^V(\mathbf{u}_{1,h}, \mathbf{v}_{0,h}, \mathbf{u}_{1,h}) = l_{0,h}^V(\mathbf{v}_{0,h}) & \forall \mathbf{v}_{0,h} \in \mathbf{X}_h, \\ b_h^V(\mathbf{u}_{0,h}, q_{0,h}) = 0 & \forall q_{0,h} \in W_{h,0}. \end{cases} \quad (47) $$

The discrete version of the system (19) is then written: find $(\mathbf{u}_{1,h}, p_{1,h}) \in \mathbf{X}_h \times W_{h,0}$ such that

$$ egin{cases} a_h^V(\mathbf{u}_{1,h}, \mathbf{v}_{1,h}) + \frac{1}{\rho} b_h^V(\mathbf{v}_{1,h}, p_{1,h}) + t_{hh}^V(\mathbf{u}_{0,h}, \mathbf{u}_{1,h}, \mathbf{v}_{1,h}) \\ + \frac{1}{2} t_{hh}^V(\mathbf{u}_{1,h}, \mathbf{u}_{0,h}, \mathbf{v}_{1,h}) - \frac{1}{2} t_{hh}^V(\mathbf{u}_{1,h}, \mathbf{v}_{1,h}, \mathbf{u}_{0,h}) = l_{1,h}^V(\mathbf{v}_{1,h}) & \forall \mathbf{v}_{1,h} \in \mathbf{X}_h, \\ b_h^V(\mathbf{u}_{1,h}, q_{1,h}) = 0 & \forall q_{1,h} \in W_{h,0}. \end{cases} \quad (48) $$

## 6.2 | Stability estimate for the fully discretized problem

The Euler scheme described in Section 5.2 is applied for the time discretization. By coupling it with the spatial scheme described in the previous subsection, the following systems are obtained, find $\mathbf{u}_{0,h}^{n+1}$ and $p_{0,h}^{n+1}$ (when $\mathbf{u}_{0,h}^0, \dots, \mathbf{u}_{0,h}^n$ and $p_{0,h}^0, \dots, p_{0,h}^n$ are known) such that
$$
\left\{
\begin{array}{ll}
\displaystyle \sum_{i=1}^{N_F} \int_{\omega_i} \frac{\mathbf{u}_{0,h}^{n+1} - \mathbf{u}_{0,h}^n}{\Delta t^n} \cdot \mathbf{v}_{0,h} + a_h^V(\mathbf{u}_{0,h}^{n+1}, \mathbf{v}_{0,h}) + \frac{1}{\rho} b_h^V(\mathbf{v}_{0,h}, P_{0,h}^{n+1}) & (49) \\
\displaystyle + t_{hh}^V(\mathbf{u}_{0,h}^n, \mathbf{u}_{0,h}^{n+1}, \mathbf{v}_{0,h}) + \frac{1}{2} t_{hh}^V(\mathbf{u}_{1,h}^n, \mathbf{u}_{1,h}^{n+1}, \mathbf{v}_{0,h}) & \\
\displaystyle - \frac{1}{2} t_{hh}^V(\mathbf{u}_{1,h}^n, \mathbf{v}_{0,h}, \mathbf{u}_{1,h}^{n+1}) = l_{0,h}^V(\mathbf{v}_{0,h}) & \forall \mathbf{v}_{0,h} \in \mathbf{X}_h, \\
b_h^V(\mathbf{u}_{0,h}^{n+1}, q_{0,h}) = 0 & \forall q_{0,h} \in W_{h,0}.
\end{array}
\right.
$$

find $\mathbf{u}_{1,h}^{n+1}$ and $p_{1,h}^{n+1}$ (when $\mathbf{u}_{1,h}^0, \dots, \mathbf{u}_{1,h}^n$ and $p_{1,h}^0, \dots, p_{1,h}^n$ are known) such that

$$
\left\{
\begin{array}{ll}
\displaystyle \sum_{i=1}^{N_F} \int_{\omega_i} \frac{\mathbf{u}_{1,h}^{n+1} - \mathbf{u}_{1,h}^n}{\Delta t^n} \cdot \mathbf{v}_{1,h} + a_h^V(\mathbf{u}_{1,h}^{n+1}, \mathbf{v}_{1,h}) + \frac{1}{\rho} b_h^V(\mathbf{v}_{1,h}, p_{1,h}^{n+1}) & (50) \\
\displaystyle + t_{hh}^V(\mathbf{u}_{0,h}^n, \mathbf{u}_{1,h}^{n+1}, \mathbf{v}_{1,h}) + \frac{1}{2} t_{hh}^V(\mathbf{u}_{1,h}^n, \mathbf{u}_{0,h}^{n+1}, \mathbf{v}_{1,h}) & \\
\displaystyle - \frac{1}{2} t_{hh}^V(\mathbf{u}_{1,h}^n, \mathbf{v}_{1,h}, \mathbf{u}_{0,h}^{n+1}) = l_{1,h}^V(\mathbf{v}_{1,h}) & \forall \mathbf{v}_{1,h} \in \mathbf{X}_h, \\
b_h^V(\mathbf{u}_{1,h}^{n+1}, q_{1,h}) = 0 & \forall q_{1,h} \in W_{h,0}.
\end{array}
\right.
$$

In the following Theorem, a stability estimate for the sum of the solutions $\mathbf{u}_{0,h}^{n+1}$ and $\mathbf{u}_{1,h}^{n+1}$ of the discrete sensitivity of the Navier-Stokes equations, Equation (49) and Equation (50) is provided. This stability estimate uses the Gronwall lemma and the Poincaré inequality.

**Theorem 1.** Let $\mathbf{f}_{0,h}$ and $\mathbf{f}_{1,h}$ be stationary. The FEV method is stable for the first-order sensitivity of the Navier-Stokes equations.

$$
\begin{align}
& \|\mathbf{u}_{0,h}^{n+1}\|^2 + \|\mathbf{u}_{1,h}^{n+1}\|^2 + \nu \int_0^{t_{n+1}} \|\nabla \mathbf{u}_{0,h}(s)\|^2 + \|\nabla \mathbf{u}_{1,h}(s)\|^2 ds \tag{51} \\
& \qquad + 2 \int_0^{t_n} t_{hh}^V(\mathbf{u}_{0,h}(s), \mathbf{u}_{0,h}(s + \Delta t^n), \mathbf{u}_{0,h}(s + \Delta t^n)) + t_{hh}^V(\mathbf{u}_{0,h}(s), \mathbf{u}_{1,h}(s + \Delta t^n), \mathbf{u}_{1,h}(s + \Delta t^n)) ds \nonumber \\
& \le \min(k_1(t^{n+1}), k_2(t^{n+1})) (\|\mathbf{f}_{0,h}\|^2 + \|\mathbf{f}_{1,h}\|^2), \nonumber
\end{align}
$$

with $k_1(t^{n+1}) = 12(e^{t^{n+1}} - 1)$, $k_2(t^{n+1}) = t^{n+1} \frac{48c^2}{\nu}$, and $c = \frac{\mathrm{Diam}(\Omega)}{\pi}$ the Poincaré constant.²⁸,²⁹

*Proof.* The stability estimate (51) is obtained by replacing $\mathbf{v}_{0,h}$ by $\mathbf{u}_{0,h}^{n+1}$ and $q_{0,h}$ by $p_{0,h}^{n+1}$ in Equation (49) and also replacing $\mathbf{v}_{1,h}$ by $\mathbf{u}_{1,h}^{n+1}$ and $q_{1,h}$ by $p_{1,h}^{n+1}$ in Equation (50)

$$
a_h^V(\mathbf{u}_{0,h}^{n+1}, \mathbf{u}_{0,h}^{n+1}) = a_h(\mathbf{u}_{0,h}^{n+1}, \mathbf{u}_{0,h}^{n+1}) = \sum_{i=1}^{N_F} \int_{\omega_i} -\nu \mathbf{u}_{0,h}^{n+1} \Delta \mathbf{u}_{0,h}^{n+1} \, dx = \nu \|\nabla \mathbf{u}_{0,h}^{n+1}\|^2, \quad (52)
$$

$$
\sum_{i=1}^{N_F} \int_{\omega_i} |\mathbf{u}_{0,h}^{n+1}|^2 dx = \|\mathbf{u}_{0,h}^{n+1}\|^2. \qquad (53)
$$

$$
a_h^V(\mathbf{u}_{1,h}^{n+1}, \mathbf{u}_{1,h}^{n+1}) = a_h(\mathbf{u}_{1,h}^{n+1}, \mathbf{u}_{1,h}^{n+1}) = \sum_{i=1}^{N_F} \int_{\omega_i} -\nu \mathbf{u}_{1,h}^{n+1} \Delta \mathbf{u}_{1,h}^{n+1} \, dx = \nu \|\nabla \mathbf{u}_{1,h}^{n+1}\|^2, \quad (54)
$$

$$
\sum_{i=1}^{N_F} \int_{\omega_i} |\mathbf{u}_{1,h}^{n+1}|^2 dx = \|\mathbf{u}_{1,h}^{n+1}\|^2. \qquad (55)
$$

According to the Cauchy-Schwarz and Young inequalities $\sum_{i=1}^{N_F} \int_{\omega_i} \mathbf{u}_{0,h}^n \cdot \mathbf{u}_{0,h}^{n+1} dx$ and $\sum_{i=1}^{N_F} \int_{\omega_i} \mathbf{u}_{1,h}^n \cdot \mathbf{u}_{1,h}^{n+1} dx$ are bounded as follows

$$
\sum_{i=1}^{N_F} \int_{\omega_i} \mathbf{u}_{0,h}^n \cdot \mathbf{u}_{0,h}^{n+1} \, dx \le \|\mathbf{u}_{0,h}^n\| \|\mathbf{u}_{0,h}^{n+1}\| \le \frac{1}{2} \|\mathbf{u}_{0,h}^n\|^2 + \frac{1}{2} \|\mathbf{u}_{0,h}^{n+1}\|^2, \quad (56)
$$
$$ \sum_{i=1}^{N_F} \int_{\omega_i} \mathbf{u}_{1,h}^n \cdot \mathbf{u}_{1,h}^{n+1} \, dx \le \|\mathbf{u}_{1,h}^n\| \|\mathbf{u}_{1,h}^{n+1}\| \le \frac{1}{2} \|\mathbf{u}_{1,h}^n\|^2 + \frac{1}{2} \|\mathbf{u}_{1,h}^{n+1}\|^2. \quad (57) $$

The linear forms $l_{0,h}^V(\mathbf{u}_{0,h}^{n+1})$ and $l_{1,h}^V(\mathbf{u}_{1,h}^{n+1})$ are continuous according to Appendix B. They are bounded according to the Young and Poincaré inequalities

$$ l_{0,h}^V(\mathbf{u}_{0,h}^{n+1}) \le \sqrt{12} \|\mathbf{f}_{0,h}\| \|\mathbf{u}_{0,h}^{n+1}\| \le 6 \|\mathbf{f}_{0,h}\|^2 + \frac{1}{2} \|\mathbf{u}_{0,h}^{n+1}\|^2, \quad (58) $$

$$ l_{1,h}^V(\mathbf{u}_{1,h}^{n+1}) \sqrt{12} \le \|\mathbf{f}_{1,h}\| \|\mathbf{u}_{1,h}^{n+1}\| \le \sqrt{12} \|\mathbf{f}_{1,h}\| \|\mathbf{u}_{1,h}^{n+1}\| \le 6 \|\mathbf{f}_{1,h}\|^2 + \frac{1}{2} \|\mathbf{u}_{1,h}^{n+1}\|^2. \quad (59) $$

By replacing the equalities Equation (52) and Equation (53) and the inequality Equation (56) in Equation (49) the following inequality holds

$$ \begin{aligned} \|\mathbf{u}_{0,h}^{n+1}\|^2 + \nu \Delta t^n \|\nabla \mathbf{u}_{1,h}^{n+1}\|^2 + \Delta t^n t_{hh}^V(\mathbf{u}_{0,h}^n, \mathbf{u}_{0,h}^{n+1}, \mathbf{u}_{0,h}^{n+1}) &+ \frac{1}{2} \Delta t^n t_{hh}^V(\mathbf{u}_{1,h}^n, \mathbf{u}_{1,h}^{n+1}, \mathbf{u}_{0,h}^{n+1}) \\ &+ \frac{1}{2} \Delta t^n t_{hh}^V(\mathbf{u}_{1,h}^n, \mathbf{u}_{0,h}^{n+1}, \mathbf{u}_{1,h}^{n+1}) \le \frac{1}{2} \|\mathbf{u}_{0,h}^n\|^2 + \frac{1}{2} \|\mathbf{u}_{0,h}^{n+1}\|^2 + 6 \Delta t^n \|\mathbf{f}_{0,h}\|^2. \end{aligned} \quad (60) $$

By replacing the equalities Equation (54) and Equation (55) and the inequality Equation (57) in Equation (50) the following inequality holds

$$ \begin{aligned} \|\mathbf{u}_{1,h}^{n+1}\|^2 + \nu \Delta t^n \|\nabla \mathbf{u}_{1,h}^{n+1}\|^2 + \Delta t^n t_{hh}^V(\mathbf{u}_{0,h}^n, \mathbf{u}_{1,h}^{n+1}, \mathbf{u}_{1,h}^{n+1}) \\ + \frac{1}{2} \Delta t^n t_{hh}^V(\mathbf{u}_{1,h}^n, \mathbf{u}_{0,h}^{n+1}, \mathbf{u}_{1,h}^{n+1}) - \frac{1}{2} \Delta t^n t_{hh}^V(\mathbf{u}_{1,h}^n, \mathbf{u}_{1,h}^{n+1}, \mathbf{u}_{0,h}^{n+1}) \le \frac{1}{2} \|\mathbf{u}_{1,h}^n\|^2 + \frac{1}{2} \|\mathbf{u}_{1,h}^{n+1}\|^2 + 6 \Delta t^n \|\mathbf{f}_{1,h}\|^2. \end{aligned} \quad (61) $$

By multiplying by 2 Equation (60) and summing over all the n, one has

$$ \begin{aligned} \|\mathbf{u}_{0,h}^{n+1}\|^2 + 2\nu \sum_{p=0}^{n} \Delta t^p \|\nabla \mathbf{u}_{0,h}^{p+1}\|^2 + 2 \sum_{p=0}^{n} \Delta t^p t_{hh}^V(\mathbf{u}_{0,h}^p, \mathbf{u}_{0,h}^{p+1}, \mathbf{u}_{0,h}^{p+1}) \\ + \sum_{p=0}^{n} \Delta t^p t_{hh}^V(\mathbf{u}_{1,h}^p, \mathbf{u}_{1,h}^{p+1}, \mathbf{u}_{0,h}^{p+1}) - \sum_{p=0}^{n} \Delta t^p t_{hh}^V(\mathbf{u}_{1,h}^p, \mathbf{u}_{0,h}^{p+1}, \mathbf{u}_{1,h}^{p+1}) \le \sum_{p=0}^{n} \Delta t^p \|\mathbf{u}_{0,h}^{n+1}\|^2 + 12 \sum_{p=0}^{n} \Delta t^p \|\mathbf{f}_{0,h}\|^2. \end{aligned} \quad (62) $$

By multiplying by 2 Equation (61) and summing all over the n, one has

$$ \begin{aligned} \|\mathbf{u}_{1,h}^{n+1}\|^2 + 2\nu \sum_{p=0}^{n} \Delta t^p \|\nabla \mathbf{u}_{1,h}^{p+1}\|^2 + 2 \sum_{p=0}^{n} \Delta t^p t_{hh}^V(\mathbf{u}_{0,h}^p, \mathbf{u}_{1,h}^{p+1}, \mathbf{u}_{1,h}^{p+1}) \\ + \sum_{p=0}^{n} \Delta t^p t_{hh}^V(\mathbf{u}_{1,h}^p, \mathbf{u}_{0,h}^{p+1}, \mathbf{u}_{1,h}^{p+1}) - \sum_{p=0}^{n} \Delta t^p t_{hh}^V(\mathbf{u}_{1,h}^p, \mathbf{u}_{1,h}^{p+1}, \mathbf{u}_{0,h}^{p+1}) \le \sum_{p=0}^{n} \Delta t^p \|\mathbf{u}_{1,h}^{p+1}\|^2 + 12 \sum_{p=0}^{n} \Delta t^p \|\mathbf{f}_{1,h}\|^2. \end{aligned} \quad (63) $$

By summing Equation (62) and Equation (63), one has

$$ \begin{aligned} \|\mathbf{u}_{0,h}^{n+1}\|^2 + \|\mathbf{u}_{1,h}^{n+1}\|^2 + 2\nu \sum_{p=0}^{n} \Delta t^p (\|\nabla \mathbf{u}_{0,h}^{p+1}\|^2 + \|\nabla \mathbf{u}_{1,h}^{p+1}\|^2) \\ + 2 \sum_{p=0}^{n} \Delta t^p t_{hh}^V(\mathbf{u}_{0,h}^p, \mathbf{u}_{0,h}^{p+1}, \mathbf{u}_{0,h}^{p+1}) + 2 \sum_{p=0}^{n} \Delta t^p t_{hh}^V(\mathbf{u}_{0,h}^p, \mathbf{u}_{1,h}^{p+1}, \mathbf{u}_{1,h}^{p+1}) \\ \le \sum_{p=0}^{n} \Delta t^p \|\mathbf{u}_{0,h}^{p+1}\|^2 + \sum_{p=0}^{n} \Delta t^p \|\mathbf{u}_{1,h}^{p+1}\|^2 + 12 \sum_{p=0}^{n} \Delta t^p (\|\mathbf{f}_{0,h}\|^2 + \|\mathbf{f}_{1,h}\|^2). \end{aligned} \quad (64) $$
Applying the rectangle method to Equation (64) gives

$$ 
\begin{aligned}
& \|\mathbf{u}_{0,h}^{n+1}\|^2 + \|\mathbf{u}_{1,h}^{n+1}\|^2 + 2\nu \int_0^{t^{n+1}} \|\nabla \mathbf{u}_{0,h}^{n+1}(s)\|^2 + \|\nabla \mathbf{u}_{1,h}^{n+1}(s)\|^2 ds \\
& + 2\int_0^{t_n} t_{hh}^V(\mathbf{u}_{0,h}(s), \mathbf{u}_{0,h}(s + \Delta t^n), \mathbf{u}_{0,h}(s + \Delta t^n)) + t_{hh}^V(\mathbf{u}_{0,h}(s), \mathbf{u}_{1,h}(s + \Delta t^n), \mathbf{u}_{1,h}(s + \Delta t^n)) ds \\
& \le \int_0^{t^{n+1}} (\|\mathbf{u}_{0,h}(s)\|^2 + \|\mathbf{u}_{1,h}(s)\|^2) ds + 12t^{n+1}(\|\mathbf{f}_{0,h}\|^2 + \|\mathbf{f}_{1,h}\|^2).
\end{aligned}
\quad (65) $$

By setting $z(t^{n+1}) = \int_0^{t^{n+1}} \|\mathbf{u}_{0,h}^{n+1}(s)\|^2 + \|\mathbf{u}_{1,h}^{n+1}(s)\|^2 ds$, then replacing it in Equation (65) and multiplying by $e^{-t^{n+1}}$, one has

$$ e^{-t^{n+1}}(z'(t^{n+1}) - z(t^{n+1})) \le 12e^{-t^{n+1}}t^{n+1}(\|\mathbf{f}_{0,h}\|^2 + \|\mathbf{f}_{1,h}\|^2). \quad (66) $$

Integrating in time Equation (66) and multiplying it by $e^{t^{n+1}}$, gives

$$ z(t^{n+1}) \le 12(e^{t^{n+1}} - 1 - t^{n+1})(\|\mathbf{f}_{0,h}\|^2 + \|\mathbf{f}_{1,h}\|^2). \quad (67) $$

Then, $z(t)$ is replaced in Equation (65), and one has

$$ 
\begin{aligned}
& \|\mathbf{u}_{0,h}^{n+1}\|^2 + \|\mathbf{u}_{1,h}^{n+1}\|^2 + \nu \int_0^{t^{n+1}} \|\nabla \mathbf{u}_{0,h}^{n+1}(s)\|^2 + \|\nabla \mathbf{u}_{1,h}^{n+1}(s)\|^2 ds \\
& + 2\int_0^{t_n} t_{hh}^V(\mathbf{u}_{0,h}(s), \mathbf{u}_{0,h}(s + \Delta t^n), \mathbf{u}_{0,h}(s + \Delta t^n)) + t_{hh}^V(\mathbf{u}_{0,h}(s), \mathbf{u}_{1,h}(s + \Delta t^n), \mathbf{u}_{1,h}(s + \Delta t^n)) ds \\
& \le k_1(t^{n+1})(\|\mathbf{f}_{0,h}\|^2 + \|\mathbf{f}_{1,h}\|^2),
\end{aligned}
\quad (68) $$

with $k_1(t^{n+1}) = 12(e^{t^{n+1}} - 1)$. The two linear forms $l_{0,h}^V(\mathbf{u}_{0,h}^{n+1})$ and $l_{1,h}^V(\mathbf{u}_{1,h}^{n+1})$ can also be bounded using the Poincaré inequality.

$$ l_{0,h}(\mathbf{u}_{0,h}^{n+1}) \le c\sqrt{12}\|\mathbf{f}_{0,h}\|\|\nabla \mathbf{u}_{0,h}^{n+1}\| \le \frac{24c^2}{\nu}\|\mathbf{f}_{0,h}\|^2 + \frac{\nu}{2}\|\nabla \mathbf{u}_{0,h}^{n+1}\|^2, \quad (69) $$

$$ l_{1,h}(\mathbf{u}_{1,h}^{n+1}) \le c\sqrt{12}\|\mathbf{f}_{1,h}\|\|\nabla \mathbf{u}_{1,h}^{n+1}\| \le \frac{24c^2}{\nu}\|\mathbf{f}_{1,h}\|^2 + \frac{\nu}{2}\|\nabla \mathbf{u}_{1,h}^{n+1}\|^2. \quad (70) $$

By plugging the inequalities (52), (56), and (69) in Equation (49) the following inequality holds

$$ 
\begin{aligned}
& \frac{1}{2}\|\mathbf{u}_{0,h}^{n+1}\|^2 + \frac{\nu}{2}\Delta t^n\|\nabla \mathbf{u}_{0,h}^{n+1}\|^2 + \Delta t^n t_{hh}^V(\mathbf{u}_{0,h}^n, \mathbf{u}_{0,h}^{n+1}, \mathbf{u}_{0,h}^{n+1}) \\
& + \frac{1}{2}\Delta t^n t_h^V(\mathbf{u}_{1,h}^n, \mathbf{u}_{1,h}^{n+1}, \mathbf{u}_{0,h}^{n+1}) - \frac{1}{2}\Delta t^n t_h^V(\mathbf{u}_{1,h}^n, \mathbf{u}_{0,h}^{n+1}, \mathbf{u}_{1,h}^{n+1}) \le \frac{1}{2}\|\mathbf{u}_{0,h}^n\|^2 + \frac{24c^2}{\nu}\Delta t^n\|\mathbf{f}_{0,h}\|^2.
\end{aligned}
\quad (71) $$

By plugging the inequalities (54), (57), and (70) in Equation (50) the following inequality holds

$$ 
\begin{aligned}
& \frac{1}{2}\|\mathbf{u}_{1,h}^{n+1}\|^2 + \frac{\nu}{2}\Delta t^n\|\nabla \mathbf{u}_{1,h}^{n+1}\|^2 + \Delta t^n t_{hh}^V(\mathbf{u}_{0,h}^n, \mathbf{u}_{1,h}^{n+1}, \mathbf{u}_{1,h}^{n+1}) \\
& + \frac{1}{2}\Delta t^n t_h^V(\mathbf{u}_{1,h}^n, \mathbf{u}_{0,h}^{n+1}, \mathbf{u}_{1,h}^{n+1}) - \frac{1}{2}\Delta t^n t_h^V(\mathbf{u}_{1,h}^n, \mathbf{u}_{1,h}^{n+1}, \mathbf{u}_{0,h}^{n+1}) \le \frac{1}{2}\|\mathbf{u}_{1,h}^n\|^2 + \frac{24c^2}{\nu}\Delta t^n\|\mathbf{f}_{1,h}\|^2.
\end{aligned}
\quad (72) $$

By multiplying by 2 Equation (71) and Equation (72), and summing over all the $n$ the two equations, then by summing this two obtained equations, we have the following inequality
$$
\begin{align}
& \|\mathbf{u}_{0,h}^{n+1}\|^2 + \|\mathbf{u}_{1,h}^{n+1}\|^2 + \nu \sum_{p=0}^{n} \Delta t^p (\|\nabla \mathbf{u}_{0,h}^{p+1}\|^2 + \|\nabla \mathbf{u}_{1,h}^{p+1}\|^2) \nonumber \\
& + 2 \sum_{p=0}^{n} \Delta t^p t_{hh}^V (\mathbf{u}_{0,h}^p, \mathbf{u}_{0,h}^{p+1}, \mathbf{u}_{0,h}^{p+1}) + 2 \sum_{p=0}^{n} \Delta t^p t_{hh}^V (\mathbf{u}_{0,h}^p, \mathbf{u}_{1,h}^{p+1}, \mathbf{u}_{1,h}^{p+1}) \tag{73} \\
& \leq \sum_{p=0}^{n} \Delta t^p \frac{48c^2}{\nu} (\|\mathbf{f}_{0,h}\|^2 + \|\mathbf{f}_{1,h}\|^2). \nonumber
\end{align}
$$

Applying the rectangle method to Equation (73) gives

$$
\begin{align}
& \|\mathbf{u}_{0,h}^{n+1}\|^2 + \|\mathbf{u}_{1,h}^{n+1}\|^2 + \nu \int_0^{t^{n+1}} \|\nabla \mathbf{u}_{0,h}^{n+1}(s)\|^2 + \|\nabla \mathbf{u}_{1,h}^{n+1}(s)\|^2 \, ds \nonumber \\
& + 2 \int_0^{t_n} t_{hh}^V (\mathbf{u}_{0,h}(s), \mathbf{u}_{0,h}(s + \Delta t^n), \mathbf{u}_{1,h}(s + \Delta t^n)) + t_{hh}^V (\mathbf{u}_{0,h}(s), \mathbf{u}_{1,h}(s + \Delta t^n), \mathbf{u}_{1,h}(s + \Delta t^n)) \, ds \tag{74} \\
& \leq k_2(t^{n+1}) (\|\mathbf{f}_{0,h}\|^2 + \|\mathbf{f}_{1,h}\|^2), \nonumber
\end{align}
$$

with $k_2(t^{n+1}) = \frac{48c^2}{\nu} t^{n+1}$. Finally, according to Equation (68) and Equation (74)

$$
\begin{align}
& \|\mathbf{u}_{0,h}^{n+1}\|^2 + \|\mathbf{u}_{1,h}^{n+1}\|^2 + \nu \int_0^{t^{n+1}} \|\nabla \mathbf{u}_{0,h}(s)\|^2 + \|\nabla \mathbf{u}_{1,h}(s)\|^2 \, ds \nonumber \\
& + 2 \int_0^{t_n} t_{hh}^V (\mathbf{u}_{0,h}(s), \mathbf{u}_{0,h}(s + \Delta t^n), \mathbf{u}_{0,h}(s + \Delta t^n)) + t_{hh}^V (\mathbf{u}_{0,h}(s), \mathbf{u}_{1,h}(s + \Delta t^n), \mathbf{u}_{1,h}(s + \Delta t^n)) \, ds \tag{75} \\
& \leq \min(k_1(t^{n+1}), k_2(t^{n+1})) (\|\mathbf{f}_{0,h}\|^2 + \|\mathbf{f}_{1,h}\|^2). \nonumber
\end{align}
$$

7 | NUMERICAL RESULTS

This section presents a numerical test already performed in Reference 11 where the sensitivity equations, the average,
and the variance of the uncertain variables are computed using the Taylor expansion method.

7.1 | Uncertainty propagation

This section aims to show how sensitivity, more precisely, the PCM, can be used for a first-order estimate of the
variance of the model output, and how the estimated variance can, in turn, be used to compute the Confidence
intervals (CI) for the physical variables of the Navier-Stokes Equation (1a)-(1d) (i.e., the horizontal and vertical
velocity and the pressure), using the PCM. The CI provides a range in which the true value of the uncertain vari-
ables may lie with a prescribed probability. In this context, the parameter *a* is a random variable with a known
distribution, expected value *μ* and standard deviation *σ*. Let *Y*(**x**, *t*; *a*) be a physical variable, (the solution of the
Navier-Stokes Equation (1a)-(1d)), with *μ*<sub>Y</sub> its expected value and *σ*<sub>Y</sub> its standard deviation. In the absence of knowl-
edge about the distribution of *Y*, the Chebyshev inequality<sup>11,33</sup> is employed and the CI of *Y* is consequently defined
as follows

$$
CI_Y = \left[ \mu_Y - \frac{\sigma_Y}{\sqrt{\alpha}}, \mu_Y + \frac{\sigma_Y}{\sqrt{\alpha}} \right], \qquad (76)
$$

where $\alpha$ denotes the confidence level.

We must first estimate the mean $\mu_Y$ and the standard deviation $\sigma_Y$ to compute such a CI. To do that, the PCM is used, and the variable $Y(\mathbf{x}, t, a)$ is written by its PCE. The mean and the standard deviation are calculated as follows
$$ 
\begin{aligned}
Y(\mathbf{x}, t; a) &= Y_0(\mathbf{x}, t)\psi_0(a) + Y_1(\mathbf{x}, t)\psi_1(a) \\
&= Y_0(\mathbf{x}, t) + \left(\frac{a - \mu}{\sigma}\right) Y_1(\mathbf{x}, t).
\end{aligned} 
\quad (77) $$

Using the expansion Equation (77), the average of Y is given

$$ \mu_Y(\mathbf{x}, t) = E(Y(\mathbf{x}, t; a)) = E\left(Y_0(\mathbf{x}, t) + \frac{a-\mu}{\sigma}Y_1(\mathbf{x}, t)\right) = Y_0(\mathbf{x}, t) = \tilde{Y}_0(\mathbf{x}, t) + \sigma \tilde{Y}_1(\mathbf{x}, t), \quad (78) $$

and the variance is obtained

$$ \sigma_Y^2(\mathbf{x}, t) = E\left((Y(\mathbf{x}, t; a) - Y_0(\mathbf{x}, t))^2\right) = Y_1^2(\mathbf{x}, t) = \sigma^2 \tilde{Y}_1^2(\mathbf{x}, t). \quad (79) $$

One can remark that the computation of the mean and the variance of $Y(\mathbf{x}, t, a)$ are performed with just two simulations, one for the state Equation (20) and another one for the sensitivity Equation (21).

## 7.2 | Flow past a square cylinder

The domain $\Omega$ is defined as in Figure 4. In this numerical test, we have an uncertainty on the inflow boundary condition $\mathbf{u} = -\frac{4a}{\ell^2}y(\ell - y)\mathbf{n}$. We consider $a$ to be normally distributed with average $\mu = 0.25$ and standard deviation $\sigma = 7.5 \times 10^{-3}$. Equation (20) and (21) are solved with homogeneous Dirichlet boundary conditions on $\Gamma_{top} \cup \Gamma_{bottom} \cup \Gamma_{obst}$, homogeneous Neumann boundary conditions on $\Gamma_{out}$, and uncertain inflow boundary condition

$$ \tilde{\mathbf{u}}_0 = -\frac{4}{\ell^2}(\mu + \sigma)y(\ell - y)\mathbf{n} \quad \text{and} \quad \tilde{\mathbf{u}}_1 = -\frac{4}{\ell^2}y(\ell - y)\mathbf{n} \quad \text{on } \Gamma_{in}. $$

The equations are solved on a mesh with a spatial step varying from 0.005 to 0.01. The discretization of both the state and the sensitivity equations is developed and solved in TrioCFD. The results obtained by the PCM are compared with those obtained using Taylor expansion and Monte Carlo methods.

In Figures 5–7, we compare the mean estimated using the PCE, the Monte Carlo method and the Taylor expansion and we show them on a horizontal cross-section ($y = 0.2$) and a vertical one ($x = 0.6$). A similarly thorough procedure is followed for the standard deviation, and the curves are presented in Figures 6–8. As one can see, the three methods give very similar results.

In Figures 9 and 10, we show the confidence intervals computed according to (76) for $\alpha = 0.05$. The first-order polynomial chaos method and the Monte Carlo one give very similar results. The PCM is more efficient in terms of computational time (2 simulations) if compared the Monte Carlo approach, which requires 1300 simulations (for the Monte Carlo estimates to reach convergence of the variance $^{11}$). The results calculated using PCM are almost identical to those calculated using Taylor expansion in Reference 11 in terms of mean and variance in this numerical test case. However, the sensitivity of the Navier–Stokes equations Equation (21) calculated using PCM are not identical to those calculated using Taylor

<figure><img src="image_8.png" /><figcaption>FIGURE 4 A rectangular channel of length $L = 2$ and width $\ell = 0.7$, with walls on the top and the bottom, and an obstacle of square section of width $d = 0.1$ at distance $x_D = 0.4$ from the inflow boundary. [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>
<figure><img src="image_9.png" /></figure>

FIGURE 5 Comparison between the mean of $u^x$, $u^y$, and $p$ calculated by the PCM, Taylor expansion, and the MC method on the horizontal cross section $y = 0.2$. [Colour figure can be viewed at wileyonlinelibrary.com]

<figure><img src="image_10.png" /></figure>

FIGURE 6 Comparison between the standard deviation of $u^x$, $u^y$, and $p$ calculated by the PCM, Taylor expansion, and the MC method on the horizontal cross section $y = 0.2$. [Colour figure can be viewed at wileyonlinelibrary.com]

<figure><img src="image_11.png" /></figure>

FIGURE 7 Comparison between the mean of $u^x$, $u^y$, and $p$ calculated by the PCM, Taylor expansion, and the MC method on the vertical cross section $x = 0.6$. [Colour figure can be viewed at wileyonlinelibrary.com]
<figure><img src="image_12.png" /></figure>

FIGURE 8 Comparison between the standard deviation of $u^x$, $u^y$, and $p$ calculated by the PCM, Taylor expansion, and the MC method on the vertical cross section $x = 0.6$. [Colour figure can be viewed at wileyonlinelibrary.com]

<figure><img src="image_13.png" /></figure>

FIGURE 9 Comparison between the confidence intervals of $u^x$, $u^y$ and $p$ computed by the PCM (blue), and the MC method (hatched red), on the horizontal cross section $y = 0.2$. [Colour figure can be viewed at wileyonlinelibrary.com]

<figure><img src="image_14.png" /></figure>

FIGURE 10 Comparison between the confidence intervals of $u^x$, $u^y$ and $p$ computed by the PCM (blue), and the MC method (hatched red), on the vertical cross section $x = 0.6$. [Colour figure can be viewed at wileyonlinelibrary.com]
expansion, since in Equation (21), there is an additional term, $2\sigma\tilde{u}_1 \cdot \nabla\tilde{u}_1$. This term is small, because it is multiplied by $\sigma$, and this explains why the numerical results are almost identical.

## 8 | CONCLUSIONS

This article proposes the PCM to model the input uncertainty and its propagation in the Navier-Stokes equations. First, the dependent variables in the state equations (1) are expressed by their PCE, and then the polynomials are calculated when the uncertain parameter $a$ is normally distributed. Second, a Galerkin projection is applied to the Navier-Stokes equations to obtain the first-order sensitivity of the Navier-Stokes equations. The sensitivity equations obtained (14) and (17) are coupled. Thus, we decoupled these equations in order to solve them numerically in a more efficient way. A discretization method called finite element-volume is also introduced and described in detail in this article. It is applied to the Navier-Stokes equations, Equation (1), and the first-order Navier-Stokes sensitivity equations Equation (21). The trilinear term was notably complex and required significant effort to be discretized accurately. A stability estimate for continuous, Equation (1) and discrete Navier-Stokes equations, Equation (43) was established. Some challenges are encountered while attempting to establish the stability of the sensitivity Navier-Stokes equations Equation (21). Nonetheless, the stability estimate of the two coupled systems Equation (14) and Equation (17) is established; the weak formulations of the coupled systems, Equation (14) and Equation (17) are discretized considering the antisymmetry of the continuous trilinear term $t^V(.,,.)$. Then, an Euler discretization in time is applied. Finally, the stability estimate is proved using the Gronwall lemma and Poincaré inequality. In conclusion, the PCM provides a reliable framework for many uncertainty quantification problems in computational fluid dynamics.

As a next step, we intend to study the convergence of the discrete sensitivity Navier-Stokes in FEV and adapt the PCM to study turbulent flows using Reynolds-Averaged Navier-Stokes (RANS) equations. Currently, we are working to address some of these challenges.

## ACKNOWLEDGMENTS

The authors would like to thank Thierry Horsin for his helpful suggestions to prove the calculation presented in Appendix B.

## DATA AVAILABILITY STATEMENT

Data available on request from the corresponding author.

## ORCID

N. Nouaime https://orcid.org/0009-0005-2649-7451

## REFERENCES

1. Bonnaire P, Pettersson P, Silva C. Intrusive generalized polynomial chaos with asynchronous time integration for the solution of the unsteady Navier-Stokes equations. *Comput Fluids*. 2021;223:104952.

2. Crestaux T, Le Maître O, Martinez J-M. Polynomial chaos expansion for sensitivity analysis. *Reliab Eng Syst Saf*. 2009;94:1161-1172.

3. Knio OM, Maître OPL. Uncertainty propagation in CFD using polynomial chaos decomposition. *Fluid Dyn Res*. 2006;38:616.

4. Knio OM, Najm HN, Ghanem RG, et al. A stochastic projection method for fluid flow: I. Basic formulation. *J Comput Phys*. 2001;173:481-511.

5. Xia L, Yuan S, Zou Z, Zou L. Uncertainty quantification of hydrodynamic forces on the DTC model in shallow water waves using CFD and non-intrusive polynomial chaos method. *Ocean Eng*. 2020;198:106920.

6. Wiener N. The homogeneous chaos. *Am J Math*. 1938;60:897-936.

7. Xiu D. *Generalized (Wiener-Askey) Polynomial Chaos*. Brown University; 2004.

8. Xiu D, Karniadakis GE. Modeling uncertainty in flow simulations via generalized polynomial chaos. *J Comput Phys*. 2003;187:137-167.

9. Mishra AA, Iaccarino G. Uncertainty estimation for Reynolds-averaged Navier-Stokes predictions of high-speed aircraft nozzle jets. *AIAA J*. 2017;55:3999-4004.

10. Vos J, Rizzi A, Darracq D, Hirschel E. Navier-Stokes solvers in European aircraft design. *Prog Aerosp Sci*. 2002;38:601-697.

11. Fiorini C, Després B, Puscas MA. Sensitivity equation method for the Navier-Stokes equations applied to uncertainty propagation. *Int J Numer Methods Fluids*. 2021;93:71-92.

12. Chalons C, Duvigneau R, Fiorini C. Sensitivity analysis for the Euler equations in Lagrangian coordinates. *International Conference on Finite Volumes for Complex Applications*. Springer; 2017:71-79.
13. Fiorini C, Chalons C, Duvigneau R. A modified sensitivity equation method for Euler equations in presence of shocks. *Numer Methods Partial Differ Equ.* 2020;36:839-867.

14. Hristova H, Etienne S, Pelletier D, Borggaard J. A continuous sensitivity equation method for time-dependent incompressible laminar flows. *Int J Numer Methods Fluids.* 2006;50:817-844.

15. Fiorini C, Puscas MA, Després B. Sensitivity analysis for a thermohydrodynamic model: uncertainty analysis and parameter estimation. *Euro J Mech B/Fluids.* 2024;105:25-33.

16. Emonot P. *Méthode de volumes éléments finis: applications aux équations de Navier-Stokes et résultats de convergence.* PhD thesis, Université Claude Bernard; 1992.

17. Fortin T. *Une méthode éléments finis à décomposition L² d'ordre élevé motivée par la simulation d'écoulement diphasique bas Mach.* PhD thesis, Université Pierre et Marie Curie; 2006.

18. Heib S. *Nouvelles discrétisations non structurées pour des écoulements de fluides à incompressibilité renforcée.* PhD thesis, Université; 2003.

19. Crouzeix M, Raviart P-A. Conforming and nonconforming finite element methods for solving the stationary stokes equations I. *ESAIM Math Modell Numer Anal.* 1973;7:33-75.

20. Angeli PE, Bieder U, Fauchet G. Overview of the TrioCFD code: Main features, V&V procedures and typical applications to nuclear engineering. *Proceedings of 16th International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH-16).* USA; 2015:-252.

21. Angeli PE, Puscas MA, Fauchet G, Cartalade A. FVCA8 benchmark for the Stokes and Navier-Stokes equations with the TrioCFD code – benchmark session. *Finite Volumes for Complex Applications VIII - Methods and Theoretical Aspects.* Springer; 2017:181-202.

22. Ciarlet PG. *Linear and Nonlinear Functional Analysis with Applications.* Vol 130. Siam Publications Library; 2013.

23. Girault V, Raviart P-A. *Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms.* Vol 5. Springer Science & Business Media; 2012.

24. Robert CP, Casella G, Casella G. *Monte Carlo statistical methods.* Vol 2. Springer; 1999.

25. Kleijnen JP. Sensitivity analysis and optimization of system dynamics models: regression analysis and statistical design of experiments. *Syst Dyn Rev.* 1995;11:275-288.

26. Paudel A, Gupta S, Thapa M, Mulani SB, Walters RW. Higher-order Taylor series expansion for uncertainty quantification with efficient local sensitivity. *Aerosp Sci Technol.* 2022;126:107574.

27. Bijl H, Lucor D, Mishra S, Schwab C. *Uncertainty Quantification in Computational Fluid Dynamics.* Lecture Notes in Computational Science and Engineering. Springer International Publishing; 2013.

28. Bebendorf M. A note on the poincaré inequality for convex domains. *Z Anal Anwend.* 2003;22:751-756.

29. Payne L, Weinberger H. An optimal poincaré inequality for convex domains. *Arch Ration Mech Anal.* 1960;5:286-292.

30. Bieder U, Graffard E. Qualification of the CFD code trio_u for full scale reactor applications. *Nucl Eng Des.* 2008;238:671-679.

31. Courant R, Friedrichs K, Lewy H. On the partial difference equations of mathematical physics. *IBM J Res Dev.* 1967;11:215-234.

32. Temam R. *Navier-Stokes equations.* by rogger temam. North-Holland, 1977. 500 pp. *J Fluid Mech.* 96(1980):827-829.

33. Jacod J, Protter P. *Probability Essentials.* Springer Science & Business Media; 2004.

**How to cite this article:** Nouaime N, Després B, Puscas MA, Fiorini C. Stability of a continuous/discrete sensitivity model for the Navier-Stokes equations. *Int J Numer Meth Fluids.* 2024;1-27. doi: 10.1002/fld.5324

## APPENDIX A. STABILITY OF DISCRETE NAVIER-STOKES EQUATIONS

In this section, the initial steps of the proof of Proposition 4 are presented. This demonstration is similar to that performed by Temam³² for the FE scheme, as the Navier-Stokes equations discretized using the FE method are the same as those discretized using the FEV method, except for the trilinear term and the linear term.

The stability estimate, Equation (44), is obtained by replacing $v_h = u_h^{n+1}$ and $q_h = p_h^{n+1}$ in Equation (43)

$$ \sum_{i=1}^{N_F} \int_{\omega_i} |\mathbf{u}_h^{n+1}|^2 d\mathbf{x} = \|\mathbf{u}_h^{n+1}\|^2, \quad a_h(\mathbf{u}_h^{n+1}, \mathbf{u}_h^{n+1}) = v \|\nabla \mathbf{u}_h^{n+1}\|^2. \qquad (A1) $$

According to the Cauchy-Schwarz and Young inequalities

$$ \sum_{i=1}^{N_F} \int_{\omega_i} \mathbf{u}_h^n \cdot \mathbf{u}_h^{n+1} d\mathbf{x} \le \|\mathbf{u}_h^n\| \|\mathbf{u}_h^{n+1}\| \le \frac{1}{2} \|\mathbf{u}_h^n\|^2 + \frac{1}{2} \|\mathbf{u}_h^{n+1}\|^2. \qquad (A2) $$
One has the continuity of the linear form, according to Appendix B,

$$l_h^V(\mathbf{u}_h^{n+1}) \le \sqrt{12} \|\mathbf{f}_h\| \|\mathbf{u}_h^{n+1}\|. \quad (A3)$$

Then according to the Young inequality, one has

$$l_h^V(\mathbf{u}_h^{n+1}) \le \sqrt{12} \|\mathbf{f}_h\| \|\mathbf{u}_h^{n+1}\| \le 6 \|\mathbf{f}_h\|^2 + \frac{1}{2} \|\mathbf{u}_h^{n+1}\|^2. \quad (A4)$$

According to the Equation (43), the following inequalities holds

$$\frac{1}{2} \|\mathbf{u}_h^{n+1}\|^2 + 2\nu \|\nabla \mathbf{u}_h^{n+1}\|^2 + \Delta t^n \|t_{hh}^V(\mathbf{u}_h^n, \mathbf{u}_h^{n+1}, \mathbf{u}_h^{n+1})\| \le \frac{1}{2} \|\mathbf{u}_h^n\|^2 + 6\Delta t^n \|\mathbf{f}_h\|^2 + \frac{1}{2} \Delta t^n \|\mathbf{u}_h^{n+1}\|^2. \quad (A5)$$

By summing Equation (A5) over all n and multiplying by 2, one has

$$\|\mathbf{u}_h^{n+1}\|^2 + 2\nu \sum_{p=0}^{n} \Delta t^p \|\nabla \mathbf{u}_h^{p+1}\|^2 + 2 \sum_{p=0}^{n} \Delta t^p t_h^V(\mathbf{u}_h^p, \mathbf{u}_h^{p+1}, \mathbf{u}_h^{p+1}) \le 12 \sum_{p=0}^{n} \Delta t^p \|\mathbf{f}_h\|^2 + \sum_{p=0}^{n} \Delta t^p \|\mathbf{u}_h^{p+1}\|^2. \quad (A6)$$

Applying the rectangle method to the Equation (A6), the following inequality holds

$$\begin{aligned} & \|\mathbf{u}_h^{n+1}\|^2 + 2\nu \int_0^{t^{n+1}} \|\nabla \mathbf{u}_h(s)\|^2 ds + 2 \int_0^{t^n} t_{hh}^V(\mathbf{u}_h(s), \mathbf{u}_h(s + \Delta t), \mathbf{u}_h(s + \Delta t)) ds \\ & \le 12t^{n+1} \|\mathbf{f}_h\|^2 + \int_0^{t^{n+1}} \|\mathbf{u}_h(s)\|^2 ds. \end{aligned} \quad (A7)$$

Setting $z(t^{n+1}) = \int_0^{t^{n+1}} \|\mathbf{u}_h(s)\|^2 ds$ then replacing in Equation (A7) and multiplying by $e^{-t^{n+1}}$ gives

$$(e^{-t^{n+1}} z(t^{n+1}))' \le 12e^{-t^{n+1}} t^{n+1} \|\mathbf{f}_h\|^2. \quad (A8)$$

Integrating in time Equation (A8) and multiplying it by $e^{t^{n+1}}$, one has

$$z(t^{n+1}) \le 12 \|\mathbf{f}_h\|^2 (e^{t^{n+1}} - 1 - t^{n+1}). \quad (A9)$$

Finally, $z(t^{n+1})$ is replaced in Equation (A7), and the following inequality holds

$$\begin{aligned} & \|\mathbf{u}_h^{n+1}\|^2 + 2\nu \int_0^{t^{n+1}} \|\nabla \mathbf{u}_h(s)\|^2 ds + 2 \int_0^{t^n} t_{hh}^V(\mathbf{u}_h(s), \mathbf{u}_h(s + \Delta t), \mathbf{u}_h(s + \Delta t)) ds \le k_1(t^{n+1}) \|\mathbf{f}_h\|^2, \\ & \|\mathbf{u}_h^{n+1}\|^2 + \nu \int_0^{t^{n+1}} \|\nabla \mathbf{u}_h(s)\|^2 ds + 2 \int_0^{t^n} t_{hh}^V(\mathbf{u}_h(s), \mathbf{u}_h(s + \Delta t), \mathbf{u}_h(s + \Delta t)) ds \le k_1(t^{n+1}) \|\mathbf{f}_h\|^2, \end{aligned} \quad (A10)$$

with $k_1(t^{n+1}) = 12(e^{t^{n+1}} - 1)$. The linear form $l_h^V(\mathbf{u}_h^{n+1})$ can also be bounded according to Cauchy-Schwarz and Poincaré inequalities as follows

$$l_h^V(\mathbf{u}_h^{n+1}) \le \sqrt{12} \|\mathbf{f}_h\| \|\mathbf{u}_h^{n+1}\| \le c \sqrt{12} \|\mathbf{f}_h\| \|\nabla \mathbf{u}_h^{n+1}\| \le \frac{24c^2}{\nu} \|\mathbf{f}_h\|^2 + \frac{\nu}{2} \|\nabla \mathbf{u}_h^{n+1}\|^2. \quad (A11)$$

By plugging the inequalities (A1), (A2), and (A11) in Equation (43), the following inequalities holds

$$\frac{1}{2} \|\mathbf{u}_h^{n+1}\|^2 + \frac{\nu}{2} \|\nabla \mathbf{u}_h^{n+1}\|^2 + \Delta t^n \|t_{hh}^V(\mathbf{u}_h^n, \mathbf{u}_h^{n+1}, \mathbf{u}_h^{n+1})\| \le \frac{1}{2} \|\mathbf{u}_h^n\|^2 + \frac{24c^2}{\nu} \|\mathbf{f}_h\|^2. \quad (A12)$$
By multiplying Equation (A12) by 2 and summing over all the $n$, then applying the rectangle method, the following inequality is obtained

$$ \|\mathbf{u}_h^{n+1}\|^2 + \nu \int_0^{t^{n+1}} \|\nabla \mathbf{u}_h(s)\|^2 ds + 2 \int_0^{t^n} t_{hh}^V(\mathbf{u}_h(s), \mathbf{u}_h(s+\Delta t), \mathbf{u}_h(s+\Delta t)) ds \leq k_2(t^{n+1}) \|\mathbf{f}_h\|^2, \quad (A13) $$

with $k_2(t^{n+1}) = \frac{48c^2}{\nu} t^{n+1}$.

Finally, according to Equation (A10) and Equation (A13), one has

$$ \begin{aligned} & \|\mathbf{u}_h^{n+1}\|^2 + \nu \int_0^{t^{n+1}} \|\nabla \mathbf{u}_h(s)\|^2 ds + 2 \int_0^{t^n} t_{hh}^V(\mathbf{u}_h(s), \mathbf{u}_h(s+\Delta t), \mathbf{u}_h(s+\Delta t)) ds \\ & \leq \min(k_1(t^{n+1}), k_2(t^{n+1})) \|\mathbf{f}_h\|^2. \end{aligned} $$

## APPENDIX B. CONTINUITY OF THE DISCRETE LINEAR FORM (A3)

Let $A(x_1, y_1)$, $B(x_2, y_2)$, and $C(x_3, y_3)$ be three points, $G$ be the barycenter of the triangle $ABC$ and $\delta$ its area. Let $M(x, y)$ be a point of the triangle $ABC$. Then $(x, y) = \lambda(x_1, y_1) + \beta(x_2, y_2) + \gamma(x_3, y_3)$, with $\lambda, \beta, \gamma \in [0, 1]$ and $\lambda + \beta + \gamma = 1$.

Let $\mathbf{u}$ be an affine function therefore $\mathbf{u}(x, y) = \lambda\mathbf{u}(x_1, y_1) + \beta\mathbf{u}(x_2, y_2) + \gamma\mathbf{u}(x_3, y_3) = \lambda\mathbf{u}_1 + \beta\mathbf{u}_2 + \gamma\mathbf{u}_3$

$$ \begin{aligned} \int_{ABC} \mathbf{u}^2(x, y) \, dxdy &= \int_{ABC} (\lambda \mathbf{u}_1 + \beta \mathbf{u}_2 + \gamma \mathbf{u}_3)^2 \, dxdy \\ &= \int_{ABC} \lambda^2 \mathbf{u}_1^2 + \beta^2 \mathbf{u}_2^2 + \gamma^2 \mathbf{u}_3^2 + 2\lambda\beta\mathbf{u}_1\mathbf{u}_2 + 2\beta\gamma\mathbf{u}_2\mathbf{u}_3 + 2\lambda\beta\mathbf{u}_1\mathbf{u}_2 \, dxdy \\ &= \mathbf{u}_1^2 \int_{ABC} \lambda^2 \, dxdy + \mathbf{u}_2^2 \int_{ABC} \beta^2 \, dxdy + \mathbf{u}_3^2 \int_{ABC} \gamma^2 \, dxdy + 2\mathbf{u}_1\mathbf{u}_2 \int_{ABC} \lambda\beta \, dxdy \\ &\quad + 2\mathbf{u}_2\mathbf{u}_3 \int_{ABC} \beta\gamma \, dxdy + 2\mathbf{u}_1\mathbf{u}_3 \int_{ABC} \lambda\gamma \, dxdy \\ &= \delta \frac{\mathbf{u}_1^2 + \mathbf{u}_2^2 + \mathbf{u}_3^2 + \mathbf{u}_1\mathbf{u}_2 + \mathbf{u}_1\mathbf{u}_3 + \mathbf{u}_2\mathbf{u}_3}{6} = \frac{1}{2} \frac{\delta}{6} [\mathbf{u}_1^2 + \mathbf{u}_2^2 + \mathbf{u}_3^2 + (\mathbf{u}_1 + \mathbf{u}_2 + \mathbf{u}_3)^2] \\ &\geq \frac{\delta}{6} \left[ \frac{1}{2} \mathbf{u}_1^2 + \frac{1}{2} \mathbf{u}_2^2 + \frac{1}{2} \mathbf{u}_3^2 \right] \geq \frac{\delta}{6} \times \frac{1}{2} \max(|\mathbf{u}_i|^2) \quad \text{with } i \in \{1, 2, 3\}. \end{aligned} $$

Therefore $\|\mathbf{u}\|_{L^\infty(ABC)} \leq \frac{\sqrt{12}}{\sqrt{|\delta|}} \|\mathbf{u}\|_{L^2(ABC)}$. Let $\delta_i^+$ and $\delta_i^-$ be respectively the areas of the triangle $\kappa_i^+$ and $\kappa_i^-$ such that $\omega_i = \kappa_i^+ \cup \kappa_i^-$.

$$ l_h^V(\mathbf{u}_h) = \sum_{i=1}^{N_F} \left( \mathbf{u}_h(\mathbf{x}_i) \cdot \int_{\omega_i} \mathbf{f}_h \, d\mathbf{x} \right) = \sum_{i=1}^{N_F} \mathbf{u}_h(\mathbf{x}_i) \cdot \left( \int_{\kappa_i^+} \mathbf{f}_h \, d\mathbf{x} + \int_{\kappa_i^-} \mathbf{f}_h \, d\mathbf{x} \right). $$

According to the Cauchy-Schwarz inequality, one obtains

$$ \begin{aligned} l_h^V(\mathbf{u}_h) &\leq \sum_{i=1}^{N_F} \left( \frac{\sqrt{12}}{\sqrt{|\delta_i^+|}} \|\mathbf{u}_h\|_{L^2(\kappa_i^+)} \sqrt{|\delta_i^+|} \|\mathbf{f}_h\|_{L^2(\kappa_i^+)} + \frac{\sqrt{12}}{\sqrt{|\delta_i^-|}} \|\mathbf{u}_h\|_{L^2(\kappa_i^-)} \sqrt{|\delta_i^-|} \|\mathbf{f}_h\|_{L^2(\kappa_i^-)} \right) \\ &\leq \sum_{i=1}^{N_F} \sqrt{12} \|\mathbf{u}_h\|_{L^2(\omega_i)} \|\mathbf{f}_h\|_{L^2(\omega_i)} \leq \sqrt{12} \|\mathbf{u}_h\|_{L^2(\Omega)} \|\mathbf{f}_h\|_{L^2(\Omega)}, \end{aligned} $$

which is the claim (A3).