# Time-dependent semidiscrete analysis of the viscoelastic fluid flow problem using a variational multiscale stabilized formulation

Time-dependent semidiscrete analysis of the viscoelastic fluid flow problem using a variational... Abstract In this article we analyse a stabilized finite element formulation recently proposed to approximate viscoelastic fluid flows. The formulation has shown to have accuracy and robustness in the different benchmarks tested in the viscoelastic framework, and permitting the use of equal interpolation of the unknown fields. We first present results about a linearized subproblem, for which well-posedness and stability results can be proved. Then, the semidiscrete nonlinear time-dependent case is addressed using a fixed point theorem, which allows us to prove existence of a semidiscrete solution, along with error estimates. 1. Introduction Viscoelastic fluids are a specific type of nonNewtonian fluids. They are characterized by having complex and high molecular weight molecules with many internal degrees of freedom (Bird et al., 2002). The classical examples of this type of fluids are the polymer solutions and molten polymers. The basic feature of polymeric fluids is the presence of long-chain molecules. In a flow, these chain molecules are stretched out by the drag forces exerted on them by the surrounding fluid (Renardy, 2000). The natural tendency of the molecule to retract from this stretched configuration generates an elastic force, which contributes to the macroscopic stress tensor, and for this reason they are called viscoelastic fluids. The interest for fluids of this kind has increased in recent years, due to the connections with the industrial applications. This motivates the numerical and mathematical analysis of the governing equations (see, e.g., Fernández-Cara et al., 2002). For viscoelastic fluid flows, in contrast to the Navier–Stokes equations, well-posedness for general models is not well understood. For initial value problems, the existence of solutions has been proved only locally in time. Global existence in time of solutions has been proved only if the initial conditions are small perturbations of the rest state, and for the steady state case existence of solutions can be proved only for small perturbations of the Newtonian case (see Renardy, 2000 for a comprehensive review). The existence of slow steady flows of viscoelastic fluids using differential constitutive equations was proved in the study by Renardy (1985) for Hilbert spaces. In this work, the author used an iterative method to show that the solution can be bounded in a certain norm, and then he proved that all iterates converge in a weaker norm. For the time-dependent case, existence of solutions locally in time, and for small data globally in time, has been proved for Hilbert spaces in the study by Guillopé & Saut (1990). The extension to Banach spaces and a complete review of uniqueness, regularity, well-posedness and stability results can be found in Fernández-Cara et al. (2002). The existence of global weak solutions for general initial conditions using a corotational Oldroyd-B model has been proved in the study by Lions & Masmoudi (2000) using a simplification (without physical justification), which consists in replacing the velocity gradient in the stress equation by its skew-symmetric part. In the study by Barrett & Boyaval (2011) the authors proved global existence of weak solutions in two dimensions to the Oldroyd-B model regularized with the introduction of a diffusion term in the constitutive equation, and assuming homogeneous natural boundary conditions associated to this term. The proof in this paper is based on a mixed finite element interpolation of the problem. The introduction of the stress diffusion can be physically justified. An analysis of the effects it has on the numerical approximation can be found in the study by Sureshkumar & Beris (1995). From a finite element perspective, the finite element approximation of the flow of viscoelastic fluids presents several numerical difficulties. On the one hand, there are all the problems inherited from the incompressible Navier–Stokes equations, mainly the compatibility between the velocity–pressure approximation and the treatment of the nonlinear advective term. But, on top of that, now the constitutive equation is highly nonlinear, with an advective term that may lead to both global and local oscillations in the numerical approximation. Moreover, even in the case of smooth solutions, it is necessary to meet some additional compatibility conditions between the velocity and the stress interpolation in order to control velocity gradients. Elements that satisfy the compatibility requirements velocity–pressure and stress–velocity are scarce, particularly in the three-dimensional case (see, e.g., Marchal & Crochet, 1987; Fortin & Fortin, 1989; Bogaerds et al., 1999). In the study by Baaijens (1998) one can find a good review of mixed methods that satisfy the two required compatibility conditions. In the finite element framework, the work of Baranger & Sandri (1992) was one of the first where the existence of an approximate solution and error bounds were given for an Oldroyd-B fluid, using Brouwer’s fixed-point theorem, discontinuous interpolation for the elastic stresses and the method of Lesaint–Raviart for the convection of the extra stress tensor in the stationary case. Later, Sandri extended (Sandri, 1994) the study to continuous approximation of the stress field, using $$P_{1}$$-$$P_{2}$$-$$P_{1}$$ interpolation (linear-quadratic-linear) for $$\sigma$$ (stress), u (velocity) and p (pressure), respectively, and the Streamline-Upwind/Petrov-Galerkin (SUPG) method to treat the convective term in the constitutive equation. The time-dependent case of the same continuous interpolation was analysed in the study by Baranger & Wardi (1995). In the same finite element context, Picasso & Rappaz (2001) analysed a stationary nonconvective Oldroyd-B problem proving a priori and a posteriori error estimates. In this work, the authors used the Galerkin Least Squares method to stabilize the momentum equation and the Elastic Viscous Split Stress (EVSS) scheme for the constitutive equation. The extension to the time-dependent case was treated in the study by Bonito et al. (2007) for the same simplified Oldroyd-B problem, proving global existence in time in Banach spaces for small data. For a Stokes/Oldroyd-B linearized problem, Bonito & Burman (2007) presented optimal a priori error estimates using the Interior Penalty method. A similar problem was studied in the study by Ervin et al. (2005) for the steady state case, but using the Johnson–Segalman linearized constitutive model, proving existence and uniqueness of the continuous solution and of a finite element approximation under a small data assumption. In another work, Ervin & Miles (2003) analysed the Oldroyd-B time-dependent case, both in the semidiscrete and in the fully discrete cases using the SUPG method, and proving existence of a solution and deriving a priori error estimates for the numerical approximation, assuming a Taylor–Hood pair approximation for the velocity and pressure and a continuous approximation for the viscoelastic stresses. The same authors extended later the analysis to a two-fluid flow problem in the study by Ervin & Miles (2005), giving a priori error estimates for the approximation in terms of the mesh and time discretiation parameters. In the study by Pani & Yuan (2005) the authors analysed the time behaviour of the viscoelastic Oldroyd model in two dimensions using a Galerkin formulation in space; in this work, the stress is eliminated through a proper projection operator, resulting in an integro-differential equation in terms of velocity and pressure. The stabilized finite element formulation analysed in this work has its roots in the Variational Multiscale (VMS) method introduced in the study by Hughes et al. (1998) for the scalar convection–diffusion-reaction problem, and extended later to the Stokes problem in the study by Codina (2000), where the space of the subgrid scales is taken orthogonal to the finite element space. As we shall see, this is an important ingredient in the method analysed, which consists in a sort of orthogonal term-by-term stabilized formulation. The key idea behind a VMS method consists in splitting the unknowns of the problem in two scales, the finite element one and the unresolvable one, called subgrid scale. The latter needs to be approximated in a simple manner in terms of the former, so as to capture its main effect and yield a globally stable formulation for the finite element unknown, keeping therefore the number of degrees of freedom of the Galerkin method. The objective of this paper is to analyse numerically the stabilized finite element formulation proposed in the study by Castillo & Codina (2014b) for the time-dependent viscoelastic flow problem. This formulation has shown to have very good accuracy and robustness in stationary (Castillo & Codina, 2014b) and time-dependent (Castillo & Codina, 2015; Castillo et al., 2015) cases. The numerical analysis of a linearized stationary case was performed in the study by Castillo & Codina (2017b), where in addition jump functions were added to permit arbitrary discontinuous interpolations of pressure and stresses. As it is usual in the analysis of numerical approximations to flow problems, even in the Newtonian case, our analysis is based on some stringent regularity assumptions on the continuous solution. Apart from analysing a nonstandard stabilized formulation, the novelty of this work is also the treatment of some of the terms that appear in the analysis. The work is organized as follows. Section 2 defines some notation and presents general results used in the subsequent analysis. In Section 3 we present the problem to be solved and its Galerkin finite element approximation, explaining the sources of the numerical instability. Section 4 contains a description of the stabilized finite element formulation analysed. Section 5 is devoted to the numerical analysis of a linearized time-dependent subproblem, where the stability of the method and the existence and uniqueness of the solution in the semidiscrete linearized case are proved. Section 6 analyses the nonlinear case, where existence of a solution is proved using a fixed point theorem. Finally, conclusions and some remarks are summarized in Section 7. 2. Notation and preliminaries 2.1 Notation Let us introduce some notation used hereafter. As usual, given a domain $$\omega$$ of $$\mathbb{R}^{d}$$, d = 2, 3, $$W^{m,p}(\omega )$$ denotes the Sobolev space of functions whose distributional derivatives of order up to $$m\in \mathbb{N}$$ belong to $$L^{p}(\omega )$$, $$p\geqslant 1$$, endowed with the standard norm. For p = 2 we write $$W^{m,2}(\omega )\equiv H^{m}(\omega )$$. The space $${H_{0}^{1}}\left (\omega \right )$$ consists of functions in $$H^{1}\left ({\omega} \right )$$ vanishing on $$\partial \omega$$. The topological dual of $${H_{0}^{1}}\left (\omega \right )$$ is denoted by $$H^{-1}\left (\omega \right )$$, the corresponding duality pairing by $$\left \langle \cdot ,\cdot \right \rangle _{\omega }$$, and the $$L^{2}$$ inner product in $$\omega$$ (for scalars, vectors and tensors) is denoted by $$\left (\cdot ,\cdot \right )_{\omega }$$. The symbol $$\left \langle \cdot ,\cdot \right \rangle _{\omega }$$ is also used for the integral over $$\omega$$ of the product of two functions, whenever it makes sense. When $$\omega =\varOmega$$, the domain where the problem is posed, the subscript or the domain information are omitted. Referring to the norms used in the subsequent analysis, $$\left \Vert \cdot \right \Vert$$ represents the $$L^{2}\left (\varOmega \right )$$-norm, $$\left \Vert \cdot \right \Vert_{L^{p}}$$ the $$L^{p}\left (\varOmega \right )$$-norm ($$2< p <\infty$$), $$\left \Vert \cdot \right \Vert_{\infty }$$ the $$L^{\infty }\left (\varOmega \right )$$-norm, $$\left \Vert \cdot \right \Vert_{k}$$ the $$H^{k}\left (\varOmega \right )$$-norm, $$\left \Vert \cdot \right \Vert_{l,p}$$ the norm for the space $$W^{l,p}\left (\varOmega \right )$$ (in particular, $$\left \Vert \cdot \right \Vert_{1,\infty }$$ is the $$W^{1,\infty }\left (\varOmega \right )$$-norm) and $$\left \vert \cdot \right \vert_{l,p}$$ the seminorm. The symbol $$\lesssim$$ will be used for $$\leqslant$$ up to constants independent of the physical parameters and the parameters of the numerical discretization. Given a vector field v, its symmetrical gradient will be denoted by $$\nabla ^{s}v$$; it is defined as $$\nabla^{s}v=\frac{1}{2}\left(\nabla v+\left(\nabla v\right)^{T}\right)\!.$$ Finally, the temporal derivative will be written as $$\partial _{t}$$. 2.2 Preliminaries Let us introduce some basic results. The gradient of a vector field v can be bounded in terms of its symmetrical gradient as in Ciarlet (2013): \begin{align} \left\Vert \nabla v\right\Vert \leqslant\sqrt{2}\left\Vert \nabla^{s} v\right\Vert,\quad v\in\left({H_{0}^{1}}\left(\varOmega\right)\right)^{d}\!.\nonumber \end{align} Poincaré–Friedrich’s and Korn’s inequalities read as: for $$v\in \left ({H_{0}^{1}}\left (\varOmega \right )\right )^{d}$$, there exists constants $$c_{\mathrm{PF}}$$ and $$c_{\textrm{K}}$$ such that \begin{align} \left\Vert v\right\Vert \leqslant c_{\mathrm{PF}}\left\Vert \nabla v\right\Vert,\;\;\;\left\Vert v\right\Vert{}_{1} \leqslant c_{\mathrm{K}}\left\Vert \nabla^{s} v\right\Vert\!.\nonumber \end{align} From the Sobolev embedding theorems (see Girault & Raviart, 1986, for example), we can bound the $$L^{4}(\varOmega )$$-norm in terms of the $$H^{1}(\varOmega )$$-norm as follows: \begin{align} \left\Vert v\right\Vert{}_{L^{4}}\leqslant C_{04}^{12} \left\Vert v\right\Vert{}_{1}. \end{align} (2.1) Constants $$c_{\mathrm{PF}}$$, $$c_{\mathrm{K}}$$ and $$C_{04}^{12}$$ depend on the shape and size of the domain $$\varOmega$$. We will use $${{c}_{\mathrm{PF}}}$$ as the constant to bound the $$L^{2}(\varOmega )$$-norm of a function by its $$H^{1}(\varOmega )$$-norm. For the finite element formulation that we shall consider, the discrete velocity field will not be pointwise divergence free and the use of the skew-symmetric counterpart of the convective term simplifies the subsequent analysis. Then, for u, v, $${w}\in ({H^{1}_{0}}(\varOmega ))^{d}$$ we define $$\tilde{c}\left({w},u, v\right)= \frac{1}{2}\left(c\left({w},u, v\right)-c\left({w}, v,u\right)\right),\quad \textrm{with}\quad c\left({w},u, v\right)= \langle{w}\cdot\nabla u, v\rangle.$$ The following properties are satisfied by this skew-symmetric form: $$\tilde{c}\left ({w},u, v\right )=c\left ({w},u, v\right )$$ when ∇⋅ w = 0 and either w = 0 or u = 0 or v = 0 on $$\partial \varOmega$$. $$\tilde{c}({w},u,u)=0$$, $${w},u\in (H^{1}(\varOmega ) )^{d}$$. $$\tilde{c}({w},u, v) \leqslant \left (C_{04}^{12}\right )^{2} \Vert{w}\Vert _{1} \Vert u\Vert _{1}\Vert v\Vert _{1}$$, $${w},u, v\in (H^{1}(\varOmega ))^{d}$$. $$\tilde{c}\left ({w},u, v\right ) \leqslant \Vert{w}\Vert _{\infty } ( \Vert u\Vert _{1}\Vert v\Vert + \Vert u\Vert \Vert v\Vert _{1})$$, $${w}\in (L^{\infty }(\varOmega ))^{d}$$, $$u, v\in (H^{1}(\varOmega ))^{d}$$. $$\tilde{c}\left ({w},u, v\right ) \leqslant \Vert v\Vert ( \Vert u\Vert _{\infty }\Vert \nabla{w}\Vert + \Vert \nabla u\Vert _{\infty }\Vert{w}\Vert )$$, $${w}\in ({H_{0}^{1}}(\varOmega ))^{d}$$, $$u\in (W^{1,\infty }(\varOmega ))^{d}$$, $$v\in (L^{2}(\varOmega ))^{d}$$. $$\tilde{c}\left ({w},u, v\right ) \leqslant \Vert v\Vert ( \Vert u\Vert \Vert \nabla{w}\Vert _{\infty } + \Vert \nabla u\Vert \Vert{w}\Vert _{\infty })$$, $${w}\in (W^{1,\infty }(\varOmega ))^{d}$$, $$u\in ({H^{1}_{0}}(\varOmega ))^{d}$$, $$v\in (L^{2}(\varOmega ))^{d}$$. Integration by parts has to be used to prove the last two properties. The same definition of the trilinear forms c and $$\tilde{c}$$ and with the same properties can be introduced when the last two arguments are tensor valued and with the appropriate regularity. Let us consider a finite element partition $$\mathscr{T}_{h} = \{ K\}$$ of the computational domain $$\varOmega$$. The diameter of an element domain $$K\in \mathscr{T}_{h}$$ is denoted by $$h_{K}$$ and the diameter of the partition, or mesh size, is defined as $$h=\max \left \{ h_{K}\mid K\in \mathscr{T}_{h}\right \}$$. We will consider for simplicity quasi-uniform families of meshes, and thus all the element diameters can be bounded above and below by constants multiplying h. Defining $$V_{h}:=\left \{ v_{h} : \varOmega \rightarrow \mathbb{R}~\vert ~ v_{h}\vert _{K}\in P_{k}\left (K\right )\:\forall \ K\in \mathscr{T}_{h}\right \}$$, $$P_{k} (K)$$ being the set of polynomials of degree k in K, we can write the following inverse inequality: there is a constant $$c_{\mathrm{inv}}$$, independent of the mesh size h, such that \begin{align} \left\Vert v_{h}\right\Vert{}_{l,\,p,\,K} \leqslant c_{\mathrm{inv}}h^{m-l+\min\left(0,\frac{d}{p}-\frac{d}{q}\right)}\left\Vert v_{h}\right\Vert{}_{m,\,q,\,K}, \end{align} (2.2) for all $$l\geqslant m, l\in [1,+\infty ]$$, and all finite element functions $$v_{h}$$ defined on $$K\in \mathscr{T}_{h}$$ (see Ern & Guermond, 2004, for example), where $$\left \Vert v_{h}\right \Vert_{l,\,p,\,K}$$ is the norm of $$v_{h}$$ in $$W^{l,\,p}(K)$$. For $$v\in W^{l,\,p}\left (\varOmega \right )$$ we have the following interpolation estimate: there exists a constant C independent of h, such that for all $$0\leqslant l\leqslant k+1$$, $$1\leqslant p\leqslant \infty$$, there holds \begin{align} \left\Vert v-P\left[v\right]\right\Vert{}_{L\,^{p}}\leqslant Ch^{l}\left| v \right|{}_{l,\,p}, \end{align} (2.3) where $$P\left [ v \right ]\in V_{h}$$ is the $$L^{2}\left (\varOmega \right )$$-projection of v in $$V_{h}$$ (see Ern & Guermond, 2004 for more details). From Sobolev’s embedding and the stability of P on finite element spaces we also have that \begin{align} \left\Vert \nabla\left(v-P\left[v\right]\right)\right\Vert{}_{\infty} \leqslant C_{\infty,3}\left\Vert v \right\Vert_{3}\!, \end{align} (2.4) where $$C_{\infty ,3}$$ is a positive constant independent of h. 3. Problem statement and Galerkin finite element discretization 3.1 The boundary value problem Let $$\varOmega$$ be an open set of $$\mathbb{R}^{d}$$ occupied by the fluid, assumed to be bounded and polyhedral, and let $$\partial \varOmega$$ be its boundary. Additionally, consider the time interval $$\left ]0,T\right [$$, with $$T<\infty$$. The incompressible and isothermal viscoelastic fluid flow problem can be written as $$\rho\left({\partial_{t}u}+u\cdot\nabla u\right)-\nabla\cdot\left(2\beta\mu\nabla^{s}u+\sigma\right)+\nabla p =f\qquad\textrm{in}\:\varOmega,\:t\in\left]0,T\right[,$$ (3.1) \begin{align} \nabla\cdot u &=0\qquad\textrm{in}\:\varOmega,\:t\in\left]0,T\right[, \end{align} (3.2) \begin{align} \frac{1}{2\mu}\sigma-(1-\beta)\nabla^{s} u+\frac{\lambda}{2\mu}\left({\partial_{t}\sigma}+ u\cdot\nabla\sigma\right)\nonumber\end{align} (3.3) \begin{align} -\frac{\lambda}{2\mu}\left(\sigma\cdot\nabla u+\left(\nabla u\right)^{T}\cdot\sigma\right) &={0}\qquad\textrm{in}\:\varOmega,\:t\in\left]0,T\right[,\end{align} (3.3) \begin{align} u &={0}\qquad\textrm{on}\:\partial\varOmega,\:t\in\left]0,T\right[,\end{align} (3.4) \begin{align} u\mid_{t=0}\ = u_{0}\qquad\!\!\textrm{in}\:\varOmega, \end{align} (3.5) \begin{align} \sigma\mid_{t=0}\ =\sigma_{0}\qquad\!\!\textrm{in}\:\varOmega. \end{align} (3.6) The unknowns of the problem are the velocity field $$u\left ({x},t\right )$$, the pressure $$p\left ({x},t\right )$$ and the viscoleastic or elastic part $$\sigma \left ({x},t\right )$$ of the extra stress tensor. The physical parameters are the dynamic viscosity $$\mu$$, the density of the fluid $$\rho$$, a real parameter $$\beta \in \left [0,1\right ]$$ to define the amount of viscous or solvent viscosity $$\left (\mu _{s}=\beta \mu \right )$$ and elastic or polymeric viscosity $$\left (\mu _{p}=(1-\beta )\mu \right )$$ in the fluid, and the relaxation time $$\lambda$$ that represents the elasticity of the fluid. Finally, $$f\in (H^{-1}(\varOmega ))^{d}$$ is the external volume force applied to the fluid confined in $$\varOmega$$. For viscoelastic fluids, the problem is incomplete without the definition of a constitutive equation for the elastic stresses $$\left (\sigma \right )$$. A large variety of approaches exist to define it (see Bird et al., 1987a,b for a complete description). In this work, we use the classical Oldroyd-B constitutive model (3.3) for this purpose. The conservation laws (3.1)–(3.2) together with the Oldroyd-B constitutive equation (3.3) are a mixed parabolic-hyperbolic problem that needs to be complemented with initial (3.5)–(3.6) and boundary (3.4) conditions to close the problem. For simplicity in the exposition, we will consider the simplest boundary condition for the velocity and no boundary conditions for the stress field. With respect to the initial conditions, $$u_{0}$$ and $$\sigma _{0}$$ are functions defined on the whole domain $$\varOmega$$. For a complete description of the mathematical structure of the problem, we refer to Fernández-Cara et al. (2002) and Renardy (1989). 3.2 Variational form To write the weak form of problem (3.1)–(3.3) we need to introduce some functional spaces. Let $${\mathscr{V}}=({H_{0}^{1}}(\varOmega ))^{d}$$, $${\varUpsilon }:=\left \{ \tau \mid \tau \in (L^{2}\left (\varOmega \right ))_{\mathrm{sym}}^{d\times d},~{w}\cdot \nabla \tau \in (L^{2}(\varOmega ))^{d\times d}_{\mathrm{sym}} ~\forall \,{w}\in{\mathscr{V}}\right \}$$ (subscript sym standing for symmetric tensors) and $$Q=L^{2}(\varOmega )/\mathbb{R}$$, the spaces of the velocity, the elastic stresses and the pressure, respectively. If we denote $$U:=\left (u,p,\sigma \right )$$, $${\mathscr{X}}:={\mathscr{V}}\times Q\times{\varUpsilon }$$, the weak form of the problem consists in finding $$U:\left ]0,T\right [\rightarrow{\mathscr{X}}$$ such that the initial conditions are satisfied and \begin{align} \rho\left({\partial_{t} u}, v\right)+\rho\langle u\cdot\nabla u, v\rangle+\left(2\beta\mu\nabla^{s} u+\sigma,\nabla^{s} v\right)-\left(p,\nabla\cdot v\right) =\left\langle f, v\right\rangle\!, \end{align} (3.7) \begin{align}\qquad\qquad\qquad\qquad\qquad \left(\nabla\cdot u,q\right) =0, \end{align} (3.8) \begin{align} \left(\frac{1}{2\mu}\sigma-(1-\beta)\nabla^{s} u,\tau\right)+\frac{\lambda}{2\mu}\left({\partial_{t}\sigma}+ u\cdot\nabla\sigma,\tau\right)-\frac{\lambda}{2\mu}\left(g\left(u,\sigma\right),\tau\right) =0, \end{align} (3.9) for all $${V}=\left (v,q,\tau \right )\in{\mathscr{X}}$$. The last term of the constitutive equation represents the traction or rotational term, defined as $$g\left (u,\sigma \right ):=\sigma \cdot \nabla u+\left (\nabla u\right )^{T}\cdot \sigma$$. In a compact form, problem (3.7)–(3.9) can be written as: \begin{align} \left(\mathscr{D}_{t}\left(U\right),{V}\right)+B\left(u,\sigma;U,{V}\right)=\left\langle f, v\right\rangle\!, \end{align} (3.10) for all $${V}\in{\mathscr{X}}$$, where \begin{align} B\left(\hat{u},\hat{\sigma};U,{V}\right)= &\, \rho\langle\hat{u}\cdot\nabla u, v\rangle+2\beta\mu\left(\nabla^{s} u,\nabla^{s} v\right)+\left(\sigma,\nabla^{s} v\right)-\left(p,\nabla\cdot v\right)+\left(\nabla\cdot u,q\right)\\ & +\frac{1}{2\mu}\left(\sigma,\tau\right)-(1-\beta)\left(\nabla^{s} u,\tau\right)+\frac{\lambda}{2\mu}\left(\hat{u}\cdot\nabla\sigma,\tau\right)-\frac{\lambda}{2\mu}\left(g\left(\hat{u},\hat{\sigma}\right),\tau\right)\!,\nonumber \end{align} (3.11) and \begin{align} \left(\mathscr{D}_{t}\left(U\right),{V}\right):=\rho\left({\partial_{t} u}, v\right)+\frac{\lambda}{2\mu}\left({\partial_{t}\sigma},\tau\right)\!. \end{align} (3.12) 3.3 Galerkin finite element discretization From $$\mathscr{T}_{h}$$ we may construct conforming finite element spaces for the velocity, the pressure and the elastic stress, $${\mathscr{V}}_{h}\subset{\mathscr{V}}$$, $$Q_{h}\subset Q$$ and $${\varUpsilon }_{h}\subset{\varUpsilon }$$, respectively. Denoting $${\mathscr{X}}_{h}={\mathscr{V}}_{h}\times Q_{h}\times{\varUpsilon }_{h}$$, the Galerkin finite element approximation of problem (3.10) consists in finding $$U_{h}:\left ]0,T\right [\rightarrow{\mathscr{X}}_{h}$$ such that \begin{align} \left(\mathscr{D}_{t}\left(U_{h}\right),{V}_{h}\right)+B\left(u_{h},\sigma_{h};U_{h},{V}_{h}\right)=\left\langle \,f, v_{h}\right\rangle\!, \end{align} (3.13) for all $${V}_{h}\in{\mathscr{X}}_{h}$$, and satisfying the appropriate initial conditions. Until now, we have posed no restrictions on the choice of the finite element spaces. However, let us analyse the numerical stability of problem (3.13). If we take $${V}_{h}=U_{h1}=\left ((1-\beta ) u_{h},(1-\beta )p_{h},\sigma _{h}\right )$$, it is found that \begin{align} B\left({u}_{h},{\sigma}_{h};U_{h},U_{h1}\right) =2(1-\beta)\beta\mu\left\Vert \nabla^{s} u_{h}\right\Vert{}^{2}+\frac{1}{2\mu}\left\Vert \sigma_{h}\right\Vert{}^{2} -\frac{\lambda}{2\mu}\left(g(u_{h},\sigma_{h}),\sigma_{h}\right)\!. \end{align} (3.14) It is seen from (3.14) that B is not coercive in $${\mathscr{X}}_{h}$$, and we can ensure only control on $$\left \Vert \sigma _{h}\right \Vert$$ for all $$\beta$$ assuming $$\lambda \nabla u_{h}$$ to be small enough. To ensure the control of $$p_{h}$$ and $$\nabla ^{s} u_{h}$$, one has then to choose finite element spaces satisfying \begin{align} \underset{q_{h}\in{{Q}}_{h}}{\textrm{inf}}\,\underset{v_{h}\in{\mathscr{V}}_{h}}{\textrm{sup}}\frac{\left(q_{h},\nabla\cdot v_{h}\right)}{\left\Vert v_{h}\right\Vert{}_{{\mathscr{V}}_{h}}\left\Vert q_{h}\right\Vert{}_{{{Q}}_{h}}} \geqslant C_{1}, \end{align} (3.15) \begin{align} \underset{v_{h}\in{\mathscr{V}}_{h}}{\textrm{inf}}\,\underset{\tau_{h}\in\mathscr{{\varUpsilon}}_{h}}{\textrm{sup}}\frac{\left(\tau_{h},\nabla^{s} v_{h}\right)}{\left\Vert \tau_{h}\right\Vert{}_{\mathscr{{\varUpsilon}}_{h}}\left\Vert v_{h}\right\Vert{}_{{\mathscr{V}}_{h}}} \geqslant C_{2}, \end{align} (3.16) where $$C_{1}$$ and $$C_{2}$$ are positive constants. It is therefore required that the finite element spaces satisfy (3.15) and (3.16), which is a stringent requirement inherited from the mixed form of the Navier–Stokes problem (Castillo & Codina, 2014a). The two compatibility conditions of the viscoelastic flow problem do not allow us the use of an arbitrary interpolation for the different fields because the scheme may become unstable. The implementation of inf–sup stable elements is a possible solution for this problem; however, from the numerical point of view, the spaces that fulfil these conditions are limited and complex, particularly when the problem needs to be solved in three dimensions. The other possibility is to use a stabilized formulation that permits the use of any interpolation for the variables, which is the approach studied in this work. Note that the constitutive equation is of convective nature, and therefore, some kind of stabilization technique has to be used even if inf–sup stable elements are used, and likewise for the momentum equation. 4. Stabilized finite element method In general, a stabilized formulation consists of replacing B in (3.10) by another multilinear form $$B_{\textrm{stab}}$$, possibly mesh dependent, designed to enhance stability without upsetting accuracy. The formulation we analyse, proposed in the study by Castillo & Codina (2014b), is described next. For the sake of conciseness, let us consider equal order continuous interpolation for all variables. The case of discontinuous pressures and stresses can be treated using the technique employed in the study by Castillo & Codina (2017b). Thus, let us consider $${\mathscr{X}}_{h}={\mathscr{V}}_{h}\times Q_{h}\times{\varUpsilon }_{h}$$, where $${\mathscr{V}}_{h}=\left [V_{h}\cap{H_{0}^{1}}(\varOmega )\right ]^{d}$$, $$Q_{h}=[V_{h}\cap \mathscr{C}^{0}(\varOmega )]/\mathbb{R}$$ and $${\varUpsilon }_{h}=\left [V_{h}\cap \mathscr{C}^{0}(\varOmega )\right ]_{\textrm{sym}}^{d\times d}$$. The method consists in replacing (3.13) by the following problem: find $$U_{h} :\ ]0, T[\ \rightarrow{\mathscr{X}}_{h}$$ such that \begin{align} \left(\mathscr{D}_{t}\left(U_{h}\right),{V}_{h}\right)+B_{\textrm{stab}}\left(u_{h},\sigma_{h};U_{h},{V}_{h}\right)=\left\langle \,f, v\right\rangle\!, \end{align} (4.1) for all $${V}_{h}\in{\mathscr{X}}_{h}$$, where \begin{align} B_{\textrm{stab}}\left(\hat{u}_{h},\hat{\sigma}_{h};U_{h},{V}_{h}\right)=B\left(\hat{u}_{h},\hat{\sigma}_{h};U_{h},{V}_{h}\right)+B^{*}\left(\hat{u}_{h},\hat{\sigma}_{h};U_{h},{V}_{h}\right)\!, \end{align} (4.2) and $$B^{*}$$represents the additional stabilization terms added to the Galerkin formulation. Using the same notation as in the study by Castillo & Codina (2014b), we can define $$B^{*}$$ as $$B^{*}\left(\hat{u}_{h},\hat{\sigma}_{h};U_{h},{V}_{h}\right)= S_{1}^{\perp}\left(\hat{u}_{h};U_{h},{V}_{h}\right)+S_{2}^{\perp}\left(U_{h},{V}_{h}\right)+S_{3}^{\perp}\left(\hat{u}_{h},\hat{\sigma}_{h};U_{h},{V}_{h}\right)\!,$$ where \begin{align} S_{1}^{\perp}\left(\hat{u}_{h};U_{h},{V}_{h}\right)= &\, \alpha_{u}(\hat{u}_{h})\left(P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla u_{h}\right],P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla v_{h}\right]\right)+\alpha_{u}(\hat{u}_{h})\left(P_{u}^{\perp}\left[\nabla p_{h}\right],P_{u}^{\perp}\left[\nabla q_{h}\right]\right)\nonumber\\ & +\alpha_{u}(\hat{u}_{h})\left(P_{u}^{\perp}\left[\nabla\cdot\sigma_{h}\right],P_{u}^{\perp}\left[(1-\beta)\nabla\cdot\tau_{h}\right]\right)\!, \end{align} (4.3) \begin{align} S_{2}^{\perp}\left(U_{h},{V}_{h}\right)= \alpha_{p}(\hat{u}_{h})\left(P_{p}^{\perp}\left[\nabla\cdot u_{h}\right],P_{p}^{\perp}\left[\nabla\cdot v_{h}\right]\right)\!, \end{align} (4.4) \begin{align} S_{3}^{\perp}\left(\hat{u}_{h},\hat{\sigma}_{h};U_{h},{V}_{h}\right)=\, & \alpha_{\sigma}(\hat{u}_{h})\left(P_{\sigma}^{\perp}\left[-(1-\beta)\nabla^{s} u_{h}+\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla\sigma_{h}-g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right)\right]\right.\!,\nonumber\\ & \left.P_{\sigma}^{\perp}\left[-\nabla^{s} v_{h}+\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla\tau_{h}+g^{*}\left(\hat{u}_{h},\tau_{h}\right)\right)\right]\right)\!, \end{align} (4.5) and $$g^{*}\left (\hat{u}_{h},\tau _{h}\right )$$ represents the adjoint operator of $$g\left (\hat{u}_{h},\tau _{h}\right )$$, defined as $$g^{*}\left (\hat{u}_{h},\tau _{h}\right ):=\tau _{h}\cdot \left (\nabla \hat{u}_{h}\right )^{T}+ \ \nabla \hat{u}_{h}\cdot \tau _{h}$$. Here $$P_{u}^{\perp }=I-P_{u}$$, where $$P_{u}:L^{2}(\varOmega )\rightarrow \tilde{{\mathscr{V}}}_{h}$$ is the $$L^{2}(\varOmega )$$-projection onto $$\tilde{{\mathscr{V}}}_{h}$$, the velocity space without boundary conditions, and $$P_{p}^{\perp }$$ and $$P_{\sigma }^{\perp }$$ are defined in an analogous way. We will also need $$P_{u,0}:=L^{2}(\varOmega )\rightarrow{{\mathscr{V}}}_{h}$$, that is to say, the $$L^{2}(\varOmega )$$-projection onto the velocity space with boundary conditions. Note that all the stabilization terms in (4.3)–(4.5) are multiplied by $$\alpha _{i}$$, $$i=u,p,\sigma$$. These terms are the components of the stabilization parameter matrix, that can be defined as $${\alpha}=\operatorname{diag}\left(\alpha_{u}{I}_{d},\alpha_{p},\alpha_{\sigma}{I}_{d\times d}\right)\!,$$ where \begin{align} \alpha_{u}(\hat{u}_{h})= \left[c_{1}\frac{\mu}{h^{2}}+c_{2}\frac{\rho\left\Vert \hat{u}_{h}\right\Vert{}_{\infty}}{h}\right]^{-1}, \end{align} (4.6) \begin{align} \alpha_{p}(\hat{u}_{h})= \frac{h^{2}}{c_{1}\alpha_{u}(\hat{u}_{h})}, \end{align} (4.7) \begin{align} \alpha_{\sigma}(\hat{u}_{h})= \left[c_{3}\frac{1}{2\mu}+\frac{\lambda}{2\mu}\left(c_{4}\frac{\left\Vert \hat{u}_{h}\right\Vert{}_{\infty}}{h}+2c_{5}\left\Vert \nabla\hat{u}_{h}\right\Vert{}_{\infty}\right)\right]^{-1}, \end{align} (4.8) where $$c_{i}$$, i = 1, 2, 3, 4, 5, are algorithmic constants and $${I}_{d}$$ and $${I}_{d\times d}$$ are the second and fourth order identity tensors, respectively. A general approach to design the terms of the stabilization parameter matrix was proposed in the study by Codina (2009) for the three-field Stokes problem. In this work, it is shown that the parameters can be uniquely determined by dimensionality, assuming that this matrix is diagonal. In general, to get optimal control the stabilization parameters need to be evaluated elementwise (or even pointwise). The constant expression adopted in this work allows us to simplify the analysis. The stabilizing mechanisms introduced by the terms $$S_{1}^{\perp }$$, $$S_{2}^{\perp }$$ and $$S_{3}^{\perp }$$ are the following. The first component of $$S_{1}^{\perp }$$ gives control on the convective term, the second component gives control on the pressure gradient and the third term gives control on the divergence of the viscoelastic stress. The term $$S_{2}^{\perp }$$ is not a must, but in some cases it improves stability of the problem. Finally, the term $$S_{3}^{\perp }$$ adds stability in the constitutive equation. Note that some of the components of this last term are the convective–convective term of the viscoelastic stress tensor and an equivalent EVSS-structure component, amongst other cross local inner product terms (see Castillo & Codina, 2014b for more details of this spatial stabilized formulation). The addition of these three terms permits the approximation of convection dominated problems both in velocity and in stress and the implementation of equal order interpolation for all the unknowns. The orthogonal projections introduce consistency errors, but of optimal order, a key point in the design of accurate nonresidual-based methods. For stationary problems, the resulting formulation turns out to have optimal order of convergence, as checked numerically in the study by Castillo & Codina (2014b) for linear and quadratic elements. The term-by-term form of $$S_{1}^{\perp }$$ was proposed instead of a residual-based one because the former shows a better numerical behaviour in problems where high gradients in pressure and stress are present (see Castillo & Codina, 2017a for more details about this fact). We will need a condition on the interpolating spaces that holds in the case of equal order interpolations (see Codina & Blasco, 2000) and that can be written as follows: given $${a}_{h}, v_{h}\in{{\mathscr{V}}}_{h}$$, $$q_{h}\in{{{Q}}}_{h}$$, $$\tau _{h}\in{\varUpsilon }_{h}$$ and $${z}_{h}:=\rho{a}_{h}\cdot \nabla v_{h}+\nabla q_{h}-\nabla \cdot \tau _{h}$$, there holds \begin{align} \left\Vert{z}_{h}\right\Vert \leqslant c_{\mathrm{m}}\left(\left\Vert P_{u,0}\left({z}_{h}\right)\right\Vert +\left\Vert P_{u}^{\perp}\left({z}_{h}\right)\right\Vert \right), \end{align} (4.9) for a constant $$c_{\mathrm{m}}>0$$. According to this condition, the component of $$P_{u}\left ({z}_{h}\right )$$ that corresponds to the boundary of $$\varOmega$$ can be bounded in terms of the right-hand side (RHS) of equation (4.9). To prove this, one can use the macroelement technique employed in the study by Codina & Blasco (2000). 5. Linearized time-dependent case The numerical analysis of the stabilized formulation presented in this work is divided in two steps. In this section we present the stability analysis of the linearized case. The second part (Section 6), is devoted to the nonlinear analysis. 5.1 Linearized stabilized semidiscrete problem The equations for incompressible viscoelastic flows have several nonlinear terms, both in the momentum and in the constitutive equation. In the former we have the convective term and in the latter we have the term corresponding to the convection of stresses and the rotational term arising from the objective derivative of stresses. On top of that, the stabilization terms depend also on the velocity, introducing therefore additional nonlinearities. As it is usual for incompressible flow problems, for the convective term of the momentum equation we can use a Picard scheme as linearization technique, taking the advection velocity as the velocity of the previous iteration. This leads only to first order convergence, but it is a robust option. The constitutive equation is rather more complex, and sometimes the implementation of combined algorithms is recommended. See for example the work of Castillo & Codina (2014b), where a Newton scheme was combined with a continuation method, or the work of Howell (2009), where different types of continuation methods were proposed. For simplicity in the numerical analysis, we will use only the fixed point scheme for all the nonlinear terms, including the terms associated to the matrix of stabilization parameters $${\alpha }$$. Therefore, we analyse the following semidiscrete linearized problem: given $$\hat{u}_{h}:\ ]0,T[\ \rightarrow{\mathscr{V}}_{h}$$ and $$\hat{\sigma }_{h}: \ ]0,T[\ \rightarrow{\varUpsilon }_{h}$$, find $$U_{h}:\left ]0,T\right [\rightarrow{\mathscr{X}}_{h}$$ such that \begin{align} \rho\left({\partial_{t} u_{h}},v_{h}\right)&+\rho\tilde{c}\left(\hat{u}_{h}, u_{h}, v_{h}\right)+\left(2\beta\mu\nabla^{s}u_{h}+\sigma_{h},\nabla^{s} v_{h}\right)-\left(p_{h},\nabla\cdot v_{h}\right)\nonumber\\ &+\alpha_{u}\left(P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla u_{h}\right],P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla v_{h}\right]\right)+\alpha_{p}\left(P_{p}^{\perp}\left[\nabla\cdot u_{h}\right],P_{p}^{\perp}\left[\nabla\cdot v_{h}\right]\right)\nonumber \\ &+\alpha_{\sigma}\left(P_{\sigma}^{\perp}\left[(1-\beta)\nabla^{s}u_{h}-\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla\sigma_{h}-g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right)\right],P_{\sigma}^{\perp}\left[\nabla^{s} v_{h}\right]\right) =\left\langle\, f, v_{h}\right\rangle\!, \end{align} (5.1) \begin{align} \qquad\qquad\qquad\qquad\qquad\qquad\;\left(\nabla\cdot u_{h},q_{h}\right)+\alpha_{u}\left(P_{u}^{\perp}\left[\nabla p_{h}\right],P_{u}^{\perp}\left[\nabla q_{h}\right]\right) =0, \end{align} (5.2) \begin{align} \left(\frac{1}{2\mu}\sigma_{h}-(1-\beta)\nabla^{s} u_{h},\tau_{h}\right)&+\frac{\lambda}{2\mu}\left({\partial_{t}\sigma_{h}},\tau_{h}\right)+\frac{\lambda}{2\mu}\tilde{c}\left(\hat{u}_{h},\sigma_{h},\tau_{h}\right)\nonumber\\ &-\frac{\lambda}{2\mu}\left(g\left(\hat{u}_{h},\hat{\sigma}_{h}\right),\tau_{h}\right)+(1-\beta)\alpha_{u}\left(P_{u}^{\perp}\left[\nabla\cdot\sigma_{h}\right],P_{u}^{\perp}\left[\nabla\cdot\tau_{h}\right]\right)\nonumber \\ &+\alpha_{\sigma}\left(P_{\sigma}^{\perp}\left[-(1-\beta)\nabla^{s}u_{h}+\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla\sigma_{h}-g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right)\right],\right.\nonumber \\ &\quad\left.P_{\sigma}^{\perp}\left[\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla\tau_{h}+g^{*}\left(\hat{u}_{h},\tau_{h}\right)\right)\right]\right) =0, \end{align} (5.3) for all $${V}_{h}=\left (v_{h},q_{h},\tau _{h}\right )\in{\mathscr{X}}_{h}$$. The initial conditions are set as appropriate projections of $$u_{0}$$ and $$\sigma _{0}$$. The stabilization parameters are computed using $$\hat{u}_{h}$$ (see (4.6)–(4.8)). 5.2 Existence and uniqueness of the semidiscrete solution The following existence and uniqueness analysis of the discrete solution was motivated by the procedure followed in the study by Burman & Fernández (2007) for the two-field Navier–Stokes problem. To prove existence and uniqueness of the discrete linearized problem (5.1)–(5.3), we shall make use of the following pressure and velocity subspaces: \begin{align*} Q_{h}^{\star}= & \left\{ q_{h}\in Q_{h}\vert\,\left(P_{u,0}^{\perp}\left[\nabla q_{h}\right],P_{u,0}^{\perp}\left[\nabla q_{h}\right]\right)=0\right\}\!,\\{\mathscr{V}}_{h}^{\operatorname{div}}= & \left\{ v_{h}\in{\mathscr{V}}_{h}\vert\,\left(q_{h},\nabla\cdot v_{h}\right)=0,\:\forall\, q_{h}\in Q_{h}^{\star}\right\}\!. \end{align*} In addition, $$Q_{h}\backslash Q_{h}^{\star }$$ will stand for the supplementary of $$Q_{h}^{\star }$$ in $$Q_{h}$$, i.e., $$Q_{h}=\left (Q_{h}\backslash Q_{h}^{\star }\right )\oplus Q_{h}^{\star }$$. To ensure that $${\mathscr{V}}_{h}^{\operatorname{div}}$$ is not trivial, we use the following lemma: Lemma 5.1 There exists a constant $$\gamma>0$$, independent of h, such that $$\underset{q_{h}\in{{Q}}_{h}^{\star}}{\textrm{inf}}\,\underset{v_{h}\in{\mathscr{V}}_{h}}{\textrm{sup}}\frac{\left(q_{h},\nabla\cdot v_{h}\right)}{\left\Vert v_{h}\right\Vert{}_{1}\left\Vert q_{h}\right\Vert}\geqslant\gamma\!.$$ Proof. Let $$q_{h}\in Q_{h}^{\star }$$. From inf–sup theory (see for example Girault & Raviart, 1986), there exists $$v_{q}\in ({H_{0}^{1}}(\varOmega ))^{d}$$ such that $$\nabla\cdot v_{q}=q_{h},\;\left\Vert v_{q}\right\Vert{}_{1}\lesssim\left\Vert q_{h}\right\Vert\!.$$ Integrating by parts, we have \begin{align*} \left\Vert q_{h}\right\Vert{}^{2}= & \left(q_{h},\nabla\cdot v_{q}\right)\\ = & -\left(\nabla q_{h}, v_{q}-P_{u,0}\left[v_{q}\right]\right)+\left(q_{h},\nabla\cdot P_{u,0}\left[v_{q}\right]\right)\\ = & \left(q_{h},\nabla\cdot P_{u,0}\left[v_{q}\right]\right)\!, \end{align*} since, as $$q_{h}\in{{Q}}_{h}^{\star }$$, then $$\nabla q_{h}\in{\mathscr{V}}_{h}$$, and $$v_{q}-P_{u,0}\left [v_{q}\right ]$$ belongs to the orthogonal to this space. In addition, using the quasi-uniformity of the mesh and the Poincaré–Friedrics inequality, we have $$\left\Vert P_{u,0}\left[v_{q}\right]\right\Vert{}_{1}^{2}\lesssim\left\Vert \nabla v_{q}\right\Vert{}^{2}\lesssim\left\Vert q_{h}\right\Vert{}^{2},$$ which completes the proof. Theorem 5.2 (Existence-uniqueness) The semidiscrete problem (5.1)–(5.3) has a unique solution. Proof. Problem (5.1)–(5.3), satisfying the initial conditions, can be written in operator form as \begin{align}{3} M_{u}\partial_{t} u_{h}+K_{u}\left(\hat{u}_{h}\right) u_{h}+Gp_{h}-D_{\sigma}\sigma_{h} =M_{u} f\; \quad \textrm{in}\:{\mathscr{V}}_{h}, \end{align} (5.4) \begin{align} D u_{h} =S_{1,p}^{\perp}p_{h}\;\; \textrm{in}\:Q_{h}, \end{align} (5.5) \begin{align} M_{\sigma}\partial_{t}\sigma_{h}+K_{\sigma}\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\sigma_{h}-S u_{h} =0\qquad\;\; \textrm{in}\:{\varUpsilon}_{h}. \end{align} (5.6) These are equations posed in terms of linear forms defined on the spaces specified in each equation. The components of these forms when applied to $$v_{h}\in{\mathscr{V}}_{h}$$, $${q}_{h}\in{{{Q}}}_{h}$$ and $$\tau _{h}\in{{\varUpsilon }}_{h}$$ are given by \begin{align*} M_{u} f (v_{h})= &\, \langle f, v_{h}\rangle,\\ (K_{u}(\hat{u}_{h}) u_{h}) (v_{h}) = &\, \rho\tilde{c}\left(\hat{u}_{h}, u_{h}, v_{h}\right)+2\beta\mu\left(\nabla^{s} u_{h}\nabla^{s} v_{h}\right)+\alpha_{u}\left(P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla u_{h}\right],P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla v_{h}\right]\right)\\ & +\alpha_{p}\left(P_{p}^{\perp}\left[\nabla\cdot u_{h}\right],P_{p}^{\perp}\left[\nabla\cdot v_{h}\right]\right)\\ & +(1-\beta)\alpha_{\sigma}\left(P_{\sigma}^{\perp}\left[(1-\beta)\nabla^{s} u_{h}-\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla\sigma_{h}-g\left(\hat{u}_{h},\sigma_{h}\right)\right)\right],P_{\sigma}^{\perp}\left[\nabla^{s} v_{h}\right]\right),\\ Gp_{h}(v_{h})= & -\left(p_{h},\nabla\cdot v_{h}\right)\!,\\ D_{\sigma}\sigma_{h}(v_{h})= & \left(\sigma_{h},\nabla^{s}v_{h}\right)\!,\\ D u_{h} (q_{h}) = & \left(\nabla\cdot u_{h},q_{h}\right)\!,\\ S_{1,p}^{\perp}p_{h} (q_{h}) = & -\alpha_{u}\left(P_{u}^{\perp}\left[\nabla q_{h}\right],P_{u}^{\perp}\left[\nabla q_{h}\right]\right)\!,\\ M_{\sigma}\partial_{t}\sigma_{h} (\tau_{h})= &\, \frac{\lambda}{2\mu}\left({\partial_{t}\sigma_{h}},\tau_{h}\right)\!,\\ (K_{\sigma}(\hat{u}_{h},\hat{\sigma}_{h})\sigma_{h})(\tau_{h})= &\, \frac{1}{2\mu}\left(\sigma_{h},\tau_{h}\right)+\frac{\lambda}{2\mu}\tilde{c}\left(\hat{u}_{h},\sigma_{h},\tau_{h}\right)+(1-\beta)\alpha_{u}\left(P_{u}^{\perp}\left[\nabla\cdot\sigma_{h}\right],P_{u}^{\perp}\left[\nabla\cdot\tau_{h}\right]\right)\\ & +\alpha_{\sigma}\left(P_{\sigma}^{\perp}\left[-(1-\beta)\nabla^{s}u_{h}+\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla\sigma_{h}-g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right)\right],\right.\\ &\quad \left.P_{\sigma}^{\perp}\left[\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla\tau_{h}+g^{*}\left(\hat{u}_{h},\tau_{h}\right)\right)\right]\right)\!,\\ S u_{h}(\tau_{h})= & -(1-\beta)\left(\nabla^{s} u_{h},\tau_{h}\right)\!. \end{align*} We also introduce the operator $$D^{1}:{\mathscr{V}}_{h} \longrightarrow ({Q}_{h}^{\star })^{\prime }$$, defined by $$D^{1}v_{h}= \left(Dv_{h}\right)_{Q_{h}^{\star}},\:\:\:\forall\, v_{h}\in{\mathscr{V}}_{h},$$ where the subscript stands for the restriction to $${Q_{h}^{\star }}$$. From Lemma 5.1, it follows that $$D^{1}$$ is surjective and $$\left (D^{1}\right )^{T}$$ is injective, and therefore, $${\mathscr{V}}_{h}^{\operatorname{div}}:=\ker \left (D^{1}\right )\neq \left \{\textbf{0}\right \}$$. Let us consider the following reduced formulation, derived from (5.1)–(5.3) with $${V}_{h}\in{\mathscr{X}}_{h}^{\star }={\mathscr{V}}_{h}^{\operatorname{div}}\times \left (Q_{h}\backslash Q_{h}^{\star }\right )\times{\varUpsilon }_{h}$$: find $$U_{h} = (u_{h},\tilde{p}_{h},\sigma _{h}) : \left ]0,T\right [\rightarrow{\mathscr{X}}_{h}^{\star }$$ such that \begin{align} M_{u}\partial_{t} u_{h}+K_{u}\left(\hat{u}_{h}\right) u_{h}+G\tilde{p}_{h}-D_{\sigma}\sigma_{h} =M_{u} f\qquad \textrm{in}\:{\mathscr{V}}_{h}^{\operatorname{div}}, \end{align} (5.7) \begin{align} D u_{h} =S_{1,p}^{\perp}\tilde{p}_{h}\quad\; \textrm{in}\:\left(Q_{h}\backslash Q_{h}^{\star}\right)\!, \end{align} (5.8) \begin{align} M_{\sigma}\partial_{t}\sigma_{h}+K_{\sigma}\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\sigma_{h}-S u_{h} =0\qquad\;\;\;\; \textrm{in}\:{\varUpsilon}_{h}. \end{align} (5.9) Since, by construction, $$Q_{h}^{\star }=\operatorname{Ker}\left (S_{1,p}^{\perp }\right )$$, we conclude that $$S_{1,p}^{\perp }$$ is invertible in $$Q_{h}\backslash Q_{h}^{\star }$$, and therefore, from (5.8) we obtain $$\tilde{p}_{h}=\left(S_{1,p}^{\perp}\right)_{Q_{h}\backslash Q_{h}^{\star}}^{-1}D u_{h}.$$ By plugging this expression into (5.7), we obtain the following equivalent problem: \begin{align} M_{u}\partial_{t} u_{h}+K_{u}\left(\hat{u}_{h}\right) u_{h}+G\left(S_{1,p}^{\perp}\right)_{Q_{h}\backslash Q_{h}^{\star}}^{-1}D u_{h}-D_{\sigma}\sigma_{h} =M_{u} f\:\:\: \textrm{in}\:{\mathscr{V}}_{h}^{\operatorname{div}}, \end{align} (5.10) \begin{align} M_{\sigma}\partial_{t}\sigma_{h}+K_{\sigma}\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\sigma_{h}-S u_{h} =0\quad\;\;\; \textrm{in}\:{\varUpsilon}_{h}, \end{align} (5.11) which is a standard Cauchy problem for $$u_{h}$$ and $$\sigma _{h}$$. Existence and uniqueness for $$u_{h}$$ and $$\sigma _{h}$$ follows by the Lipschitz continuity of $$K_{u}\left (\hat{u}_{h}\right )$$ and $$K_{\sigma }\left (\hat{u}_{h},\hat{\sigma }_{h}\right )$$. We may then recover $$\tilde{p}_{h}$$ uniquely from (5.8). Therefore, the reduced problem (5.7)–(5.9) has a unique solution. On the other hand, from (5.7) it follows that $$M_{u}\partial_{t} u_{h}+K_{u}\left(\hat{u}_{h}\right)u_{h}+G\tilde{p}_{h}-D_{\sigma}\sigma_{h}-M_{u} f \in\left(\operatorname{Ker}\left(D^{1}\right)\right)^{0},$$ with $$\left (\operatorname{Ker}\left (D^{1}\right )\right )^{0}$$ standing for the polar set of $$\operatorname{Ker}\left (D^{1}\right )$$. From Lemma 5.1, it follows that $$D^{1}$$ is an isomorphism from $$Q_{h}^{\star }$$ onto $$\left (\operatorname{Ker}\left (D^{1}\right )\right )^{0}$$. Thus, there exists a unique $$p^{1}\in Q_{h}^{\star }$$ such that \begin{align} M_{u}\partial_{t}u_{h}+K_{u}\left(\hat{u}_{h}\right) u_{h}+G\tilde{p}_{h}-D_{\sigma}\sigma_{h}-M_{u} f= Gp^{1}\:\:\:\textrm{in}\,{\mathscr{V}}_{h}. \end{align} (5.12) Therefore, from (5.12) and the reduced problem (5.7)–(5.9), it follows that the discrete problem (5.1)–(5.3) has a unique solution, given by $$\left (u_{h},\tilde{p}_{h}-p^{1},\sigma _{h}\right )$$. 5.3 Stability of the linearized semidiscrete problem We will use the following working norm in the subsequent analysis: \begin{align} \left\Vert{V}_{h}\right\Vert{}_{W}^{2}= \, & 2\beta\mu(1-\beta)\left\Vert \nabla^{s}v_{h}\right\Vert{}^{2}+\frac{1}{2\mu}\left\Vert \tau_{h}\right\Vert{}^{2}+(1-\beta)\alpha_{u}\left\Vert P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla v_{h}\right]\right\Vert^{2}\nonumber \\ & +(1-\beta)\alpha_{u}\left\Vert \rho{\partial_{t} v_{h}}+\rho\hat{u}_{h}\cdot\nabla v_{h}+\nabla q_{h}-\nabla\cdot\tau_{h}\right\Vert^{2}\nonumber\\ & +(1-\beta)\alpha_{p}\left\Vert \nabla\cdot v_{h}\right\Vert{}^{2}+(1-\beta)\alpha_{u}\left\Vert P_{u}^{\perp}\left[\nabla q_{h}\right]\right\Vert^{2}+(1-\beta)\alpha_{u}\left\Vert P_{u}^{\perp}\left[\nabla\cdot\tau_{h}\right]\right\Vert^{2}\nonumber \\ & +\alpha_{\sigma}\left\Vert \frac{\lambda}{2\mu}{\partial_{t}\tau_{h}}-(1-\beta)\nabla^{s}v_{h}+\frac{\lambda}{2\mu}\hat{u}_{h}\cdot\nabla\tau_{h}\right\Vert^{2}. \end{align} (5.13) The term multiplied by $$\alpha _{p}$$ is not strictly necessary, since it is already contained in $$\left \Vert \nabla ^{s} v_{h}\right \Vert$$, but sometimes could reinforce stability. We will keep it for generality, to see the effect of the stabilizing term associated to the pressure. From the definition of $$\Vert \cdot \Vert _{W}$$ it is observed that we have some control on the convective terms of the momentum equation and of the constitutive equation. In view of the expression of the stabilization parameters, this control remains meaningful in the convection dominated limit for the momentum equation. This property will no longer be valid in the nonlinear case, since, as it is standard in the analysis of nonlinear problems involving the Navier–Stokes equation, the results will be proved in the diffusion dominated regime. The next result states the stability of the proposed semidiscrete formulation, defined in (5.1)–(5.3). Theorem 5.3 (Stability). Let $$\left (u_{h},p_{h},\sigma _{h}\right )$$ be the solution of (5.1)–(5.3). Then, the following stability holds, for almost all $$t\in \left [0,T\right ]$$: \begin{align} & (1-\beta)\frac{\rho}{2}\left\Vert u_{h}\left(t\right)\right\Vert{}^{2}+\frac{\lambda}{4\mu}\left\Vert \sigma_{h}\left(t\right)\right\Vert{}^{2} +{\int_{0}^{t}}\left\Vert (u_{h},p_{h},\sigma_{h}) \right\Vert{}_{W}^{2}\,\mathrm{d}t \nonumber \\ & \qquad \lesssim \frac{c_{\mathrm{K}}^{2}}{2\mu}{\int_{0}^{t}}\left\Vert f\right\Vert{}_{H^{-1}}^{2}\,\mathrm{d}t+\frac{\lambda^{2}}{2\mu}{\int_{0}^{t}}\left\Vert g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right\Vert{}^{2}\,\mathrm{d}t +(1-\beta)\frac{\rho}{2}\left\Vert u_{0}\right\Vert{}^{2}+\frac{\lambda}{4\mu}\left\Vert \sigma_{0}\right\Vert{}^{2}. \end{align} (5.14) Proof. In this proof $$\varepsilon _{1},\varepsilon _{2},...$$ are positive constants used in the application of Young’s inequalities. The values will be chosen at the end of the proof. Taking $$\left (v_{h},q_{h},\tau _{h}\right )=U_{h1}=\left ((1-\beta )u_{h},(1-\beta )p_{h},\sigma _{h}\right )$$ in (5.1)–(5.3), adding up the resulting equations and using Cauchy–Schwarz’s and Young’s inequalities, we arrive at \begin{align} & (1-\beta)\frac{\rho}{2}\frac{\mathrm{d}}{\mathrm{d}t}\left\Vert u_{h}\right\Vert{}^{2}+2\mu\beta(1-\beta)\left\Vert \nabla^{s} u_{h}\right\Vert{}^{2}+\frac{\lambda}{4\mu}\frac{\mathrm{d}}{\mathrm{d}t}\left\Vert \sigma_{h}\right\Vert{}^{2}\nonumber\\ &\quad +(1-\beta)\alpha_{p}\left\Vert P_{p}^{\perp}\left[\nabla\cdot u_{h}\right]\right\Vert^{2}+(1-\beta)\alpha_{u}\left\Vert P_{u}^{\perp}\left[\nabla\cdot\sigma_{h}\right]\right\Vert^{2}\nonumber \\ & \quad+\frac{1}{2\mu}\left(1-\frac{\varepsilon_{1}}{2}-4\frac{\alpha_{\sigma}}{2\mu}\lambda^{2}\left\Vert \nabla\hat{u}_{h}\right\Vert{}_{\infty}^{2}\left(\frac{1}{2\varepsilon_{2}}+\frac{1}{2\varepsilon_{4}}\right)\right)\left\Vert \sigma_{h}\right\Vert{}^{2}\nonumber \\ & \quad+(1-\beta)\alpha_{u}\left\Vert P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla u_{h}\right]\right\Vert^{2}+(1-\beta)\alpha_{u}\left\Vert P_{u}^{\perp}\left[\nabla p_{h}\right]\right\Vert^{2}\nonumber \\ & \quad+\alpha_{\sigma}\left(1-\frac{\varepsilon_{2}}{2}-\frac{\varepsilon_{3}}{2}\right)\left\Vert P_{\sigma}^{\perp}\left[-(1-\beta)\nabla^{s} u_{h}+\frac{\lambda}{2\mu}\left({\partial_{t}\sigma_{h}}+\hat{u}_{h}\cdot\nabla\sigma_{h}\right)\right]\right\Vert^{2}\nonumber \\ &\qquad \leqslant (1-\beta)\left\langle f, u_{h}\right\rangle +\frac{\lambda^{2}}{4\mu}\left[\frac{1}{\varepsilon_{1}}\left\Vert g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right\Vert{}^{2}+\frac{\alpha_{\sigma}}{\mu}\left(\frac{1}{2\varepsilon_{3}}+\frac{\varepsilon_{4}}{2}\right)\left\Vert P_{\sigma}^{\perp}\left[g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right]\right\Vert^{2}\right]. \end{align} (5.15) If we compare the bounded terms of expression (5.15) with the terms of the working norm (5.13), we can see that the missing terms are all of them associated to the finite element space. The key point is that this missing control comes from the Galerkin part of the multilinear form $$B_{\operatorname{stab}}$$ in equation (4.2). Let us take $${V}_{h1}=\left (v_{h},q_{h},\tau _{h}\right )=\alpha _{u}\left ((1-\beta )v_{1},0,{0}\right )$$ with $$v_{1}\equiv P_{u,0}\left[\rho{\partial_{t} u_{h}}+\rho\hat{u}_{h}\cdot\nabla u_{h}+\nabla p_{h}-\nabla\cdot\sigma_{h}\right]\!.$$ Then, we can proceed in a similar way as in the $$U_{h1}$$ case. Taking $$v_{1}$$ as test function, integrating by parts the terms arising from the divergence of the total stress and the gradient of the pressure, using the fact that $$\hat{u}_{h}={0}$$ on $$\partial \varOmega$$, and applying Cauchy–Schwarz’s and Young’s inequalities and inverse inequalities, we obtain the following result \begin{align} &(1-\beta)\alpha_{u}\left(1 -\alpha_{u}c_{\operatorname{inv}}^{2}\left[\frac{\alpha_{p}}{h^{2}}\frac{\varepsilon_{6}}{2}+\frac{\rho\left\Vert \hat{u}_{h}\right\Vert{}_{\infty}}{h}\frac{\varepsilon_{7}}{2}+\frac{\alpha_{\sigma}}{h^{2}}\left(\frac{\varepsilon_{8}}{2}+\frac{\varepsilon_{9}}{2}\right)\right]\right.\nonumber\\ &\quad \left. -\beta\alpha_{u}\frac{2\mu}{h^{2}}c_{\operatorname{inv}}^{2}\frac{\varepsilon_{5}}{2}-\alpha_{u}\frac{d}{4\varepsilon_{10}}\frac{c_{\mathrm{K}}^{2}}{2\mu}\rho^{2}\left\Vert \nabla\hat{u}_{h}\right\Vert{}_{\infty}^{2}\right)\left\Vert v_{1}\right\Vert{}^{2}\nonumber \\ &\quad -2\mu(1-\beta)\left(\beta\frac{1}{2\varepsilon_{5}}+\frac{\varepsilon_{10}}{2}\right)\left\Vert \nabla^{s} u_{h}\right\Vert{}^{2}-(1-\beta)\alpha_{p}\frac{1}{2\varepsilon_{6}}\left\Vert P_{p}^{\perp}\left[\nabla\cdot u_{h}\right]\right\Vert^{2}\nonumber \\ &\quad -(1-\beta){\alpha_{u}^{2}}\frac{\rho\left\Vert \hat{u}_{h}\right\Vert{}_{\infty}}{h}\frac{1}{2\varepsilon_{7}}\left\Vert P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla u_{h}\right]\right\Vert^{2}\nonumber \\ &\quad -(1-\beta)\alpha_{\sigma}\frac{1}{2\varepsilon_{8}}\left\Vert P_{\sigma}^{\perp}\left[\frac{\lambda}{2\mu}{\partial_{t}\sigma_{h}}-(1-\beta)\nabla^{s} u_{h}+\frac{\lambda}{2\mu}\hat{u}_{h}\cdot\nabla\sigma_{h}\right]\right\Vert^{2}\nonumber \\ &\quad\quad\leqslant (1-\beta)\alpha_{u}\left\langle f, v_{1}\right\rangle +\frac{(1-\beta)}{2\varepsilon_{9}}\frac{\lambda^{2}}{2\mu}\frac{\alpha_{\sigma}}{2\mu}\left\Vert P_{\sigma}^{\perp}\left[g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right]\right\Vert^{2}. \end{align} (5.16) To obtain the control of $$P_{p}\left [\nabla \cdot u_{h}\right ]$$, we proceed taking $${V}_{h2}=(1-\beta )\alpha _{p}\left ({0},q_{2},{0}\right )$$ with $$q_{2}\equiv P_{p}\left [\nabla \cdot u_{h}\right ]$$: \begin{align} (1-\beta)\alpha_{p}\left(\left(\nabla\cdot u_{h},q_{2}\right)+\alpha_{u}\left(P_{u}^{\perp}\left[\nabla p_{h}\right],\nabla q_{2}\right)\right)\geqslant &\, (1-\beta)\alpha_{p}\left(1-c_{\operatorname{inv}}^{2}\alpha_{u}\frac{\alpha_{p}}{h^{2}}\frac{\varepsilon_{11}}{2}\right)\left\Vert P_{p}\left[\nabla\cdot u_{h}\right]\right\Vert{}^{2}\nonumber \\ & -(1-\beta)\alpha_{u}\frac{1}{2\varepsilon_{11}}\left\Vert P_{u}^{\perp}\left[\nabla p_{h}\right]\right\Vert^{2}. \end{align} (5.17) Finally, taking: $${V}_{h3}=\alpha _{\sigma }\left ({0},0,\sigma _{3}\right )$$ with $$\sigma _{3}\equiv P_{\sigma }\left [\frac{\lambda }{2\mu }{\partial _{t}\sigma _{h}}-(1-\beta )\nabla ^{s} u_{h}+\frac{\lambda }{2\mu }\hat{u}_{h}\cdot \nabla \sigma _{h}\right ]$$, and proceeding as before, we obtain \begin{align} & -\frac{1}{2\mu}\left(\frac{1}{2\varepsilon_{12}}+\frac{1}{2}\frac{1}{2\varepsilon_{17}}\right)\left\Vert \sigma_{h}\right\Vert{}^{2}-\frac{(1-\beta)}{2\varepsilon_{14}}\alpha_{u}\left\Vert P_{u}^{\perp}\left[\nabla\cdot\sigma_{h}\right]\right\Vert^{2}\nonumber\\ & +\alpha_{\sigma}\left(1-\frac{\alpha_{\sigma}}{2\mu}\frac{\varepsilon_{12}}{2}-\frac{\varepsilon_{13}}{2}-\frac{\varepsilon_{17}}{4}\frac{\alpha_{\sigma}}{2\mu}d\left(\lambda\left\Vert \nabla\hat{u}_{h}\right\Vert{}_{\infty}\right)^{2}-c_{\operatorname{inv}}^{2}(1-\beta)\alpha_{u}\frac{\alpha_{\sigma}}{h^{2}}\frac{\varepsilon_{14}}{2}\vphantom{\left[\left(\frac{\lambda\left\Vert \hat{u}_{h}\right\Vert{}_{\infty}}{2\mu h}\right)^{2}+4\left(\frac{\lambda\left\Vert \nabla\hat{u}_{h}\right\Vert{}_{\infty}}{2\mu}\right)^{2}\right]}\right.\nonumber \\ & \left.-\,c_{\operatorname{inv}}^{2}2\alpha_{\sigma}^{2}\left[\left(\frac{\lambda\left\Vert \hat{u}_{h}\right\Vert{}_{\infty}}{2\mu h}\right)^{2}+4\left(\frac{\lambda\left\Vert \nabla\hat{u}_{h}\right\Vert{}_{\infty}}{2\mu}\right)^{2}\right]\left(\frac{1}{2\varepsilon_{15}}+\frac{1}{2\varepsilon_{16}}\right)\right)\left\Vert \sigma_{3}\right\Vert{}^{2}\nonumber \\ & -\alpha_{\sigma}\frac{\varepsilon_{15}}{2}\left\Vert P_{\sigma}^{\perp}\left[\frac{\lambda}{2\mu}{\partial_{t}\sigma_{h}}-(1-\beta)\nabla^{s} u_{h}+\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla\sigma_{h}\right)\right]\right\Vert^{2}\nonumber \\ & \quad \leqslant \frac{\lambda^{2}}{4\mu}\frac{\alpha_{\sigma}}{\mu}\left(\frac{1}{2\varepsilon_{13}}\left\Vert g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right\Vert{}^{2}+\frac{\varepsilon_{16}}{2}\left\Vert P_{\sigma}^{\perp}\left[g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right]\right\Vert^{2}\right). \end{align} (5.18) Let $${V}_{h}^{*}=U_{h1}+\theta _{1}{V}_{h1}+\theta _{2}{V}_{h2}+\theta _{3}{V}_{h3}$$, with $${V}_{hi}$$, i = 1, …, 3 introduced above. Adding (5.16)–(5.18) multiplied by $$\theta _{i}$$, i = 1, …, 3, and adding also (5.15), we arrive at an expresion of the form $$LHS\left({V}_{h}^{*}\right)\leqslant RHS\left({V}_{h}^{*}\right)\!.$$ For the RHS, applying Young’s inequalities and the inverse inequality, we obtain \begin{align} RHS\left({V}_{h}^{*}\right) &\leqslant (1-\beta)\left(\frac{1}{2\mu}\left(\frac{1}{2\varepsilon_{17}}+\theta_{1}\frac{1}{2\varepsilon_{18}}\right)c_{\varOmega}^{2}\left\Vert f\right\Vert{}_{H^{-1}}^{2}\right.\nonumber\\ & \quad\left.+\ 2\mu\frac{\varepsilon_{17}}{2}\left\Vert \nabla^{s} u_{h}\right\Vert{}^{2}+\theta_{1}\alpha_{u}\left(\alpha_{u}\frac{2\mu}{h^{2}}\right)c_{\operatorname{inv}}^{2}\frac{\varepsilon_{18}}{2}\left\Vert v_{1}\right\Vert{}^{2}\right)\nonumber \\ & \quad+\frac{\lambda^{2}}{2\mu}\left(\frac{1}{2\varepsilon_{1}}+\theta_{3}\frac{\alpha_{\sigma}}{2\mu}+\frac{\alpha_{\sigma}}{2\mu}\left(\frac{1}{2\varepsilon_{3}}+\frac{\varepsilon_{4}}{2}+\theta_{1}\frac{(1-\beta)}{2\varepsilon_{9}}+\theta_{3}\frac{\varepsilon_{16}}{2}\right)\right)\left\Vert g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right\Vert{}^{2}. \end{align} (5.19) For the left-hand side (LHS), integrating inequalities (5.15), (5.16), (5.17), (5.18) and (5.19) from 0 to t, we obtain \begin{align} & (1-\beta)\frac{\rho}{2}\left\Vert u_{h}\right\Vert{}^{2}\left(t\right)+\frac{\lambda}{4\mu}\left\Vert \sigma_{h}\right\Vert{}^{2}\left(t\right)+2\mu(1-\beta)\beta C_{1}{\int_{0}^{t}}\left\Vert \nabla^{s}u_{h}\right\Vert{}^{2}\,\mathrm{d}t\nonumber\\ &\quad +\frac{1}{2\mu}C_{2}{\int_{0}^{t}}\left\Vert \sigma_{h}\right\Vert{}^{2}\,\mathrm{d}t+(1-\beta)C_{3}\alpha_{u}{\int_{0}^{t}}\left\Vert P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla u_{h}\right]\right\Vert^{2}\,\mathrm{d}t\nonumber \\ &\quad +(1-\beta)C_{4}\alpha_{p}{\int_{0}^{t}}\left\Vert P_{p}^{\perp}\left[\nabla\cdot u_{h}\right]\right\Vert^{2}\,\mathrm{d}t+(1-\beta)C_{5}\alpha_{u}{\int_{0}^{t}}\left\Vert P_{u}^{\perp}\left[\nabla\cdot\sigma_{h}\right]\right\Vert^{2}\,\mathrm{d}t\nonumber \\ &\quad +(1-\beta)C_{6}\alpha_{u}{\int_{0}^{t}}\left\Vert P_{u}^{\perp}\left[\nabla p_{h}\right]\right\Vert^{2}\,\mathrm{d}t+C_{7}\alpha_{\sigma}{\int_{0}^{t}}\left\Vert P_{\sigma}^{\perp}\left[-(1-\beta)\nabla^{s} u_{h}+\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla\sigma_{h}\right)\right]\right\Vert^{2}\,\mathrm{d}t\nonumber \\ &\quad +(1-\beta)C_{8}\alpha_{u}{\int_{0}^{t}}\left\Vert P_{u,0}\left[\rho{\partial_{t} u_{h}}+\rho\hat{u}_{h}\cdot\nabla u_{h}+\nabla p_{h}-\nabla\cdot\sigma_{h}\right]\right\Vert{}^{2}\,\mathrm{d}t\nonumber \\ &\quad +(1-\beta)C_{9}\alpha_{p}{\int_{0}^{t}}\left\Vert P_{p}\left[\nabla\cdot u_{h}\right]\right\Vert{}^{2}\,\mathrm{d}t\nonumber \\ &\quad +C_{10}\alpha_{\sigma}{\int_{0}^{t}}\left\Vert P_{\sigma}\left[\frac{\lambda}{2\mu}{\partial_{t}\sigma_{h}}-(1-\beta)\nabla^{s}u_{h}+\frac{\lambda}{2\mu}\hat{u}_{h}\cdot\nabla\sigma_{h}\right]\right\Vert^{2}\,\mathrm{d}t\nonumber \\ & \qquad \leqslant (1\!-\!\beta)\frac{c_{\varOmega}^{2}}{2\mu}C_{11}{\int_{0}^{t}}\!\left\Vert f\right\Vert{}_{H^{-1}}^{2}\,\mathrm{d}t\!+\!\frac{\lambda^{2}}{2\mu}C_{12}{\int_{0}^{t}}\!\left\Vert g\!\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\!\right\Vert{}^{2}\,\mathrm{d}t+(1-\beta)\frac{\rho}{2}\left\Vert u_{h}\left(0\right)\right\Vert{}^{2}\!+\!\frac{\lambda}{4\mu}\!\left\Vert \sigma_{h}\left(0\right)\right\Vert{}^{2}\!, \end{align} (5.20) where $$C_{i}$$, i = 1, …, 12, can be easily identified. The result then follows by choosing $$\varepsilon _{1},\ldots ,\varepsilon _{18}$$ and $$\theta _{1}, \theta _{2} ,\theta _{3}$$, in such a way that $$C_{i}>0$$ for all i, which is possible by choosing some of the constants and the algorithmic stabilization parameters sufficiently small. The missing control on $$\rho \partial _{t}u_{h}+\rho \hat{u}_{h}\cdot \nabla u_{h}+\nabla p_{h}-\nabla \cdot \sigma _{h}$$ follows from (4.9). 6. Analysis of the nonlinear problem In this section we show that under suitable conditions, a solution to the discretized system exists. Using the same procedure proposed by Ervin & Miles (2003), the proof can be subdivided in the following four steps: Define an iterative map in such a way that a fixed point of the map is a solution to the original problem. Show that the map is well defined and bounded on bounded sets. Show that there exists an invariant ball of the map. Apply Schauder’s fixed point theorem to conclude there exists a discrete solution in this ball. The results to be obtained in this section yield existence of a semidiscrete solution, as well as stability and convergence. Due to the hypotheses needed in the proof of the results (most of which are already present in the analysis of the Navier–Stokes equation), the norm in which the results are shown is weaker than the norm used in the linearized problem. In essence, we will have $$L^{\infty }(0,T)$$-control for both the $$L^{2}(\varOmega )$$-norm of velocity and stresses and $$L^{2}(0,T)$$-control for the $$H^{1}(\varOmega )$$-norm of velocity. The pressure, on the other hand, is controlled only in a norm involving the stabilization term, and not the natural $$L^{2}(\varOmega )$$-norm. The precise assumptions required on the continuous solution are: Assumption 6.1 System (3.10) has a solution $$\left (u,p,\sigma \right )$$ continuous in time and satisfying \begin{alignat}{3} & \underset{0\leqslant t\leqslant T}{\sup}\left\Vert u\right\Vert{}_{\infty}\leqslant D_{1}, \quad & \underset{0\leqslant t\leqslant T}{\sup}\left\Vert \nabla u\right\Vert{}_{\infty}\leqslant D_{2}, \quad & \underset{0\leqslant t\leqslant T}{\sup}\left\Vert u\right\Vert{}_{k+1}\leqslant D_{3},\nonumber\\ & \underset{0\leqslant t\leqslant T}{\sup}\left\Vert \sigma\right\Vert{}_{\infty}\leqslant D_{4}, \quad & \underset{0\leqslant t\leqslant T}{\sup}\left\Vert \nabla\sigma\right\Vert{}_{\infty}\leqslant D_{5}, \quad & \underset{0\leqslant t\leqslant T}{\sup}\left\Vert \sigma\right\Vert{}_{k+1}\leqslant D_{6}, \\ & \underset{0\leqslant t\leqslant T}{\sup}\left\Vert p\right\Vert{}_{k}\leqslant D_{7}, \quad & \underset{0\leqslant t\leqslant T}{\sup}\left\Vert \partial_{t} u\right\Vert{}_{k}\leqslant D_{8}, \quad & \underset{0\leqslant t\leqslant T}{\sup}\left\Vert \partial_{t}\sigma\right\Vert{}_{k}\leqslant D_{9},\nonumber \end{alignat} (6.1) for certain positive constants $$D_{i}$$i = 1, …, 9, which are supposed to be small enough. For the time-discrete problem, if $$\delta t$$ is the time step size, one usually needs a condition of the form $$\delta t \geqslant C \alpha _{u}$$ for a positive constant C, which is encountered in most stabilized finite element methods; see Codina et al. (2007) and Badia & Codina (2009) and references therein for a description of the problem and a way to avoid this restriction, which we shall not consider in this work. In the time continuous case, the boundedness in time of $$\Vert p\Vert _{k}$$ and the assumption that T is large enough allows us to prove convergence, as we show next. We do not pretend however to consider the long-term behaviour of the solution, which would require the modification of the stabilized formulation and the analysis presented in the study by Badia et al. (2010). In order to write all estimates in dimensionless form, let $$L_{d}$$ be a characteristic length of the problem and $$T_{d}$$ a characteristic time scale. These parameters may explode with the viscosity, and therefore the estimates are not valid for high Reynolds numbers. The main result on existence and convergence reads as follows: Theorem 6.2 (Convergence) Let k be the interpolation order, assumed to be the same for all variables, with $$k\geqslant d/2$$. Suppose also that Assumption 6.1 holds, that T is sufficiently large and that the $$L^{2}(0,T; (H^{-1}(\varOmega ))^{d})$$-norm of f is bounded. Then, if the viscosity is sufficiently large, there exists a solution to (5.1)–(5.3) satisfying \begin{alignat}{1} \underset{0\leqslant t\leqslant T}{\sup}\left\Vert u-u_{h}\right\Vert{}^{2}+\frac{1}{T_{d}}{\int_{0}^{T}}\left(L_{d}\left\Vert \nabla\left(u-u_{h}\right)\right\Vert \right)^{2}\,\mathrm{d}t\leqslant \, &\, u_{\star}^{2}h^{2k},\nonumber\\ \underset{0\leqslant t\leqslant T}{\sup}\left\Vert \sigma-\sigma_{h}\right\Vert{}^{2}\leqslant \, &\, \sigma_{\star}^{2}h^{2k}, \\{\int_{0}^{T}}\alpha_{u} \left\Vert \nabla\left(p-p_{h}\right)\right\Vert{}^{2}\,\mathrm{d} t\leqslant \, &\, p_{\star}^{2}h^{2k}, \nonumber \end{alignat} (6.2) where $$u_{\ast }$$, $$p_{\ast }$$ and $$\sigma _{\ast }$$ are appropriate dimensional factors that render the estimates dimensionally consistent. Proof. Step 1. The iterative map. A mapping $$\delta :L^{2}\left (0,T;{\mathscr{V}}_{h}\right )\times L^{2}\left (0,T;Q_{h}\right )\times L^{\infty }\left (0,T;{\varUpsilon }_{h}\right )\rightarrow L^{2}\left (0,T;{\mathscr{V}}_{h}\right )\times L^{2}\left (0,T;Q_{h}\right )\times L^{\infty }\left (0,T;{\varUpsilon }_{h}\right )$$ is defined via $$\left (u_{h},p_{h},\sigma _{h}\right )=\delta \left (\hat{u}_{h},\hat{p}_{h},\hat{\sigma }_{h}\right )$$, where $$U_{h} = \left (u_{h},p_{h},\sigma _{h}\right )$$ satisfies (5.1)–(5.3), for all $${V}_{h}=\left (v_{h},q_{h},\tau _{h}\right )\in{\mathscr{X}}_{h}$$. Thus, given an initial guess for the three unknowns $$\hat{U}_{h} := \left (\hat{u}_{h},\hat{p}_{h},\hat{\sigma }_{h}\right )$$, solving the above system for $$\left (u_{h},p_{h},\sigma _{h}\right )$$ gives a new approximation to the solution. Also, it is clear that the fixed point is a solution to the approximating system (5.1)–(5.3), i.e., $$\delta \left ({u}_{h},{p}_{h},{\sigma }_{h}\right )=\left ({u}_{h},{p}_{h},{\sigma }_{h}\right )$$ implies that $$\left ({u}_{h},{p}_{h},{\sigma }_{h}\right )$$ is a solution to (4.1). Step 2. The mapping $$\delta$$ is well defined and bounded on bounded sets. The existence and uniqueness of the discrete solution was proved in subsection 5.2 for a linearized problem, that can be associated to the solution of the fixed-point problem. The stability result proved in subsection 5.3 ensures that the linearized problem, that can be viewed as a fixed-point iteration of the nonlinear case, is stable and bounded under suitable regularity assumptions. Note that this step is crucial in the definition of the fixed point mapping and can be used to establish that the mapping $$\delta$$ is bounded on bounded sets. Step 3. Existence of an invariant ball. We begin defining an invariant ball. Let $$R=h^{k}$$, and for $$V = (v,q,\tau )$$ let us define the norm \begin{align} \Vert V {\Vert_{B}^{2}} := \rho \sup_{0\leqslant t\leqslant T}\Vert v\Vert^{2} + \frac{\rho}{T_{d}}{\int_{0}^{T}} (L_{d}\Vert \nabla v \Vert)^{2}\,\mathrm{d} t + \frac{\lambda}{\mu} \sup_{0\leqslant t\leqslant T} \Vert \tau\Vert^{2} + \alpha_{u} {\int_{0}^{T}} \Vert \nabla q \Vert^{2}\,\mathrm{d} t. \end{align} (6.3) As we shall see, our result could be proved in a finer norm, but the additional terms that we could include in $$\Vert \cdot \Vert _{B}$$ (see below) do not provide significant additional control for large viscosities. Let us define now the ball $$\mathscr{B}_{h}$$ as \begin{align} \mathscr{B}_{h}= \Bigl\{ V_{h} = \left(v_{h},q_{h},\tau_{h}\right) : \ ]0,T[\ \rightarrow{\mathscr{X}}_{h}~~\textrm{such that}~~ \Vert V_{h} -U\Vert_{B} \leqslant U_{\ast} R \Bigr\}, \end{align} (6.4) where $$U^{2}_{\ast } = \rho u^{2}_{\ast } + \frac{\lambda }{\mu }\sigma ^{2}_{\ast } + p^{2}_{\ast }$$ and $$u_{\ast }$$, $$p_{\ast }$$ and $$\sigma _{\ast }$$ are constructed along the proof. In the definition (6.4), $$U = \left (u,p,\sigma \right )$$ is the solution of (3.10). Let us pick now $$\hat{U}_{h} := \left (\hat{u}_{h},\hat{p}_{h},\hat{\sigma }_{h}\right )\in \mathscr{B}_{h}$$, arbitrary, and let $$U_{h} = \left (u_{h},p_{h},\sigma _{h}\right )=\delta \left (\hat{u}_{h},\hat{p}_{h},\hat{\sigma }_{h}\right )$$; it satisfies (5.1)–(5.3). Then, it is readily checked that \begin{align} \rho\left(\partial_{t}\left(u-u_{h}\right), v_{h}\right)&+\rho\left(\tilde{c}\left(u,u, v_{h}\right)-\tilde{c}\left(\hat{u}_{h},u_{h}, v_{h}\right)\right)\nonumber\\ &+2\beta\mu\left(\nabla^{s}\left(u-u_{h}\right),\nabla^{s} v_{h}\right)+\left(\sigma-\sigma_{h},\nabla^{s} v_{h}\right)-\left(p-p_{h},\nabla\cdot v_{h}\right)\nonumber \\ &-\alpha_{u}\left(P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla u_{h}\right],P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla v_{h}\right]\right)-\alpha_{p}\left(P_{p}^{\perp}\left[\nabla\cdot u_{h}\right],P_{p}^{\perp}\left[\nabla\cdot v_{h}\right]\right)\nonumber \\ &+\left(\nabla\cdot\left(u-u_{h}\right),q_{h}\right)-\alpha_{u}\left(P_{u}^{\perp}\left[\nabla p_{h}\right],P_{u}^{\perp}\left[\nabla q_{h}\right]\right)\nonumber \\ &+\frac{1}{2\mu}\left(\sigma-\sigma_{h},\tau_{h}\right)-(1-\beta)\left(\nabla^{s}\left(u-u_{h}\right),\tau_{h}\right)+\frac{\lambda}{2\mu}\left(\partial_{t}\left(\sigma-\sigma_{h}\right),\tau_{h}\right)\nonumber \\ &+\frac{\lambda}{2\mu}\left(\tilde{c}\left(u,\sigma,\tau_{h}\right)-\tilde{c}\left(\hat{u}_{h},\sigma_{h},\tau_{h}\right)\right)-\frac{\lambda}{2\mu}\left(g\left(u,\sigma\right),\tau_{h}\right)+\frac{\lambda}{2\mu}\left(g\left(\hat{u}_{h},\hat{\sigma}_{h}\right),\tau_{h}\right)\nonumber \\ &-(1-\beta)\alpha_{u}\left(P_{u}^{\perp}\left[\nabla\cdot\sigma_{h}\right],P_{u}^{\perp}\left[\nabla\cdot\tau_{h}\right]\right)\nonumber \\ &-\alpha_{\sigma}\left(P_{\sigma}^{\perp}\left[-(1-\beta)\nabla^{s}u_{h}+\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla\sigma_{h}-g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right)\right],\right.\nonumber \\ &\quad\left.P_{\sigma}^{\perp}\left[-\nabla^{s} v_{h} + \frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla\tau_{h}+g^{*}\left(\hat{u}_{h},\tau_{h}\right)\right)\right]\right) =0. \end{align} (6.5) Let us define the following approximation and interpolation errors: \begin{align*} \varLambda=u-{\mathscr{U}} & \:\:\textrm{and}\:\: E={\mathscr{U}}-u_{h},\\ \varGamma=\sigma-{\varSigma} & \:\:\textrm{and}\:\: F={\varSigma}-\sigma_{h},\\ \varPi=p-\mathscr{P} & \:\:\textrm{and}\:\:G=\mathscr{P}-p_{h}, \end{align*} where $${\mathscr{U}}=P_{u,0}\left [u\right ]$$, $${\varSigma }=P_{\sigma }\left [\sigma \right ]$$ and $$\mathscr{P}=P_{p}\left [p\right ]$$. We obviously have that $${e}_{u} = u-u_{h}=\varLambda +E$$, $${e}_{\sigma } = \sigma -\sigma _{h}=\varGamma + F$$, $$e_{p} = p-p_{h}=\varPi +G$$. Let us describe the strategy to prove that U belongs to the ball $$\mathscr{B}_{h}$$. Let $$P[U] = ({\mathscr{U}},{\varSigma }, \mathscr{P})$$. From (6.5) we will show that \begin{align} \Vert P[U] - U_{h}{\Vert_{B}^{2}} \leqslant \varphi_{1}(D) \Vert U -\hat{U}_{h}{\Vert^{2}_{B}} + \Vert U - P[U]{\Vert^{2}_{I}}, \end{align} (6.6) where $$\varphi _{1}(D)$$ is a certain function polynomial in terms of the components of the array of constants $$D = (D_{1},\dots ,D_{9})$$ introduced in Assumption 6.1 and $$\Vert \cdot \Vert _{I}$$ is a norm of what we may consider the interpolation error U − P[U]; this norm will appear along the proof. Thanks to the properties of the interpolation, we will check that \begin{align} \Vert U - P[U]{\Vert_{I}^{2}} \leqslant \varphi_{2}(D) h^{2k}, \end{align} (6.7) and we can already verify that \begin{align} \Vert U - P[U]{\Vert_{B}^{2}} \leqslant \varphi_{3}(D) h^{2k}. \nonumber \end{align} Again, $$\varphi _{2}(D)$$ and $$\varphi _{3}(D)$$ are functions polynomial in terms of the components of D. Using the triangle inequality, (6.6)–(6.7) and the fact that $$\hat{U}_{h} \in \mathscr{B}_{h}$$, we will have that \begin{align} \Vert U - U_{h}{\Vert_{B}^{2}} & \leqslant \Vert P[U] - U_{h}{\Vert_{B}^{2}} + \Vert U - P[U]{\Vert_{B}^{2}} \nonumber\\ & \leqslant \varphi_{1}(D) U^{2}_{\ast} h^{2k} + \varphi_{2}(D) h^{2k} + \varphi_{3}(D) h^{2k}.\nonumber \end{align} Taking the components of array D sufficiently small, we will be able to guarantee that \begin{align} \Vert U - U_{h}{\Vert_{B}^{2}} \leqslant U^{2}_{\ast} h^{2k},\nonumber \end{align} that is to say, $${U}_{h} \in \mathscr{B}_{h}$$. Therefore, the goal is to prove (6.6) and check that (6.7) holds. The constant $$U_{\ast }$$, from which we can obtain $$u_{\ast }$$, $$p_{\ast }$$ and $$\sigma _{\ast }$$, can be taken as $$U^{2}_{\ast } = c( \varphi _{2}(D) + \varphi _{3}(D))$$, with c > 1, provided D is such that $$\varphi _{1}(D) + c^{-1} \leqslant 1$$. Taking $$v_{h}=(1-\beta )\left (E+\gamma \alpha _{u}{P_{u,0}\left [\nabla G\right ]}\right )$$, with $$\gamma> 0$$, $$\tau _{h}= F$$ and $$q_{h}=(1-\beta )G$$ in (6.5), we obtain \begin{align} \rho(1-\beta) \frac{\mathrm{d}}{\mathrm{d}t}\left\Vert E\right\Vert{}^{2}+\frac{\lambda}{\mu} \frac{\mathrm{d}}{\mathrm{d}t}\left\Vert F\right\Vert{}^{2}+ Q(t) \lesssim N(t) + S(t) +{\overset{23}{\underset{j=1}{\sum}}R_{j}}(t), \end{align} (6.8) where \begin{align*} Q(t) & := \beta(1-\beta)\mu\left\Vert E\right\Vert{}_{1}^{2} +\frac{1}{\mu}\left\Vert F\right\Vert{}^{2}+(1-\beta)\alpha_{u}\left\Vert P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla E\right]\right\Vert^{2}+(1-\beta)\alpha_{p}\left\Vert P_{p}^{\perp}\left[\nabla\cdot E\right]\right\Vert^{2}\nonumber \\ &\qquad +(1-\beta)\alpha_{u}\left(\left\Vert P_{u}^{\perp}\left[\nabla G\right]\right\Vert{}^{2}+{\gamma\left\Vert{P_{u,0}\left[\nabla G\right]}\right\Vert{}^{2}}\right)+(1-\beta)\alpha_{u}\left\Vert P_{u}^{\perp}\left[\nabla\cdot F\right]\right\Vert^{2}\nonumber \\ & \qquad + \alpha_{\sigma}\Bigl\Vert P_{\sigma}^{\perp} \Bigl[ -(1-\beta)\nabla^{s} E+\frac{\lambda}{2\mu}\hat{u}_{h}\cdot\nabla F\Bigr] \Bigr\Vert^{2},\nonumber\\ N(t) & := -\rho(1-\beta)\left(\tilde{c}\left(u,u, E\right)-\tilde{c}\left(\hat{u}_{h},u_{h}, E\right)\right)-\frac{\lambda}{2\mu}\left(\tilde{c}\left(u,\sigma,F\right)-\tilde{c}\left(\hat{u}_{h},\sigma_{h},F\right)\right)\nonumber \\ &\qquad {-\,(1-\beta)\gamma\alpha_{u}\rho\left(\tilde{c}\left(u,u,{P_{u,0}\left[\nabla G\right]}\right)+\tilde{c}\left(\hat{u}_{h},u_{h},{P_{u,0}\left[\nabla G\right]}\right)\right)},\nonumber \\ S(t) & := -\,\alpha_{\sigma}\Bigl(P_{\sigma}^{\perp}\Bigl[-(1-\beta)\nabla^{s} E+\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla F\right)\Bigr],g^{*}(\hat{u}_{h},F)\Bigr), \end{align*} and the terms $$R_{j}(t)$$ are defined below. From (6.8), and using (4.9) for the pressure term, we have that \begin{align} \Vert P[U] - U_{h}{\Vert^{2}_{B}} & \lesssim \rho(1-\beta) \sup_{0\leqslant t\leqslant T} \Vert E\Vert^{2} + \frac{\lambda}{\mu} \sup_{0\leqslant t\leqslant T} \Vert F\Vert^{2} + {\int_{0}^{T}} Q(t) \,\mathrm{d}t \nonumber\\ &\lesssim{\int_{0}^{T}} \Biggl(N(t) + S(t) +{\overset{23}{\underset{j=1}{\sum}}R_{j}}(t) \Biggr)\,\mathrm{d}t. \end{align} (6.9) The objective is to see that all these terms in the RHS can be bounded as indicated in (6.6) or absorbed by the LHS. We will not detail the positive dimensionless constants that appear in this process, in which we will make frequent use of Young’s inequality; the parameter that appears when using it will be generically denoted by $$\varepsilon$$ or $$\xi$$, understanding that it is small enough. We could track the different appearances of these parameters as we have done in the proof of Theorem 5.3, and choose them at the end; however, since we are not interested in the values of the constants, we will proceed in a more conceptual way. We will start with the terms $$R_{j}(t)$$ of (6.8). Using the $$L^{2}$$-orthogonality between $$\varLambda$$ and E and between $$\varGamma$$ and F and Young’s inequality, we have that \begin{align*} R_{1}(t) & =-(1-\beta)\rho\left(\partial_{t}\varLambda, E\right) = 0,\\ R{}_{2}(t) & =-\frac{\lambda}{\mu}\left(\partial_{t}\varGamma,F\right) = 0,\\ R_{3}(t) & =(1-\beta)\left(\varPi,\nabla\cdot E\right)\lesssim \mu(1-\beta){\varepsilon}\left\Vert E\right\Vert{}_{1}^{2}+\frac{1}{\mu}\frac{(1-\beta)}{\varepsilon}\left\Vert \varPi\right\Vert{}^{2}. \end{align*} The term that involves the discrete error of the pressure needs a special treatment. Using (4.9) we have that \begin{align*} R_{4}(t) & =-(1-\beta)\left(\nabla\cdot\varLambda,G\right)= (1-\beta)\left(\varLambda,\nabla G\right)\\ & \lesssim \alpha_{u}^{-1}\frac{(1-\beta)}{\varepsilon}\left\Vert \varLambda\right\Vert{}^{2}+(1-\beta){\varepsilon}\alpha_{u}\left( \left\Vert{P_{u}^{\perp}\left[\nabla G\right]}\right\Vert{}^{2} + \left\Vert{{P_{u,0}\left[\nabla G\right]}}\right\Vert{}^{2} \right). \end{align*} The following terms only need Cauchy–Schwarz’s and Young’s inequalities: \begin{align*} R_{5}(t) & =-2\beta(1-\beta)\mu\left(\nabla^{s}\varLambda,\nabla^{s} E\right)\lesssim \beta(1-\beta)\mu{\varepsilon}\left\Vert E\right\Vert{}_{1}^{2}+ \beta(1-\beta)\mu\frac{1}{\varepsilon}\left\Vert \varLambda\right\Vert{}_{1}^{2},\\ R_{6}(t) & =-(1-\beta)\left(\varGamma,\nabla^{s} E\right)\lesssim \mu(1-\beta){\varepsilon}\left\Vert E\right\Vert{}_{1}^{2}+\frac{1}{\mu}(1-\beta)\frac{1}{\varepsilon}\left\Vert \varGamma\right\Vert{}^{2},\\ R_{7}(t) & =-\frac{1}{2\mu}\left(\varGamma,F\right) = 0,\\ R_{8}(t) & =(1-\beta)\left(\nabla^{s}\varLambda,F\right)\lesssim \frac{1}{\mu}(1-\beta){\varepsilon}\left\Vert F\right\Vert{}^{2} + \mu(1-\beta)\frac{1}{\varepsilon}{\left\Vert\varLambda\right\Vert{}_{1}^{2}}. \end{align*} Taking the parameter $$\varepsilon$$ small enough, it can be readily checked that the contributions from $$R_{j}$$ can be absorbed by Q(t) or are interpolation errors that behave as indicated in (6.7). Let us treat now the terms of the RHS of (6.8) that come from the stabilization. Using Cauchy–Schwarz’s and Young’s inequalities, we can easily control the following terms: \begin{align*} R_{9}(t) & = (1-\beta)\alpha_{u}\left(P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla{\mathscr{U}}\right],P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla E\right]\right),\\ R_{10}(t)\! & = (1-\beta)\alpha_{p}\left(P_{p}^{\perp}\left[\nabla\cdot{\mathscr{U}}\right],P_{p}^{\perp}\left[\nabla\cdot E\right]\right),\\ R_{11}(t) & \!= (1-\beta)\alpha_{u}\left(P_{u}^{\perp}\left[\nabla\mathscr{P}\right],P_{u}^{\perp}\left[\nabla G\right]\right),\\ R_{12}(t)\! & = (1-\beta)\alpha_{u}\left(P_{u}^{\perp}\left[\nabla\cdot{\varSigma}\right],P_{u}^{\perp}\left[\nabla\cdot F\right]\right),\\ R_{13}(t) & \!= \alpha_{\sigma}\!\left(P_{\sigma}^{\perp}\left[\!-(1-\!\beta)\nabla^{s}{\mathscr{U}}+\!\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla{\varSigma}\right)\right],P_{\sigma}^{\perp}\left[\!-(1-\beta)\nabla^{s}E+\!\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla F\!+g^{*}\left(\hat{u}_{h},F\right)\right)\!\right]\right)\!,\\ R_{14}(t) & \!= \alpha_{\sigma}\left(P_{\sigma}^{\perp}\left[-\frac{\lambda}{2\mu}g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right],P_{\sigma}^{\perp}\left[-(1-\beta)\nabla^{s} E+\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla F+g^{*}\left(\hat{u}_{h},F\right)\right)\right]\right). \end{align*} Some of these terms need additional treatment. Let us consider first $$R_{9}(t)$$. When bounding it, we need to make use of the fact that \begin{align} \alpha_{u}\left\Vert P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla{\mathscr{U}}\right]\right\Vert^{2} & \leqslant \alpha_{u}\rho^{2}\left\Vert P_{u}^{\perp}\left[\left(\hat{u}_{h}-u\right)\cdot\nabla\left({\mathscr{U}-}u\right)\right]\right\Vert^{2}+\alpha_{u}\rho^{2}\left\Vert P_{u}^{\perp}\left[u\cdot\nabla\left({\mathscr{U}-}u\right)\right]\right\Vert^{2}\nonumber\\ &\quad +\alpha_{u}\rho^{2}\left\Vert P_{u}^{\perp}\left[u\cdot\nabla u\right]\right\Vert^{2}+\alpha_{u}\rho^{2}\left\Vert P_{u}^{\perp}\left[\left(\hat{u}_{h}-u\right)\cdot\nabla u\right]\right\Vert^{2}, \end{align} (6.10) where the four terms of the RHS can be shown to have structure given by (6.6) due to Assumption 6.1 and making use of adequate interpolation estimates. Terms $$R_{10}(t)$$ to $$R_{13}(t)$$ are again bounded by terms that can be absorbed by Q(t), as well as by terms that involve the orthogonal projection of an operator applied to the projection onto the finite element space of the continuous solution. They can all be treated using the triangle inequality. For example, for the term bounding $$R_{10}(t),$$ we have \begin{align} \left\Vert P_{p}^{\perp}\left[\nabla\cdot{\mathscr{U}}\right]\right\Vert^{2} \leqslant \left\Vert P_{p}^{\perp}\left[\nabla\cdot({\mathscr{U}} - u)\right]\right\Vert^{2} + \left\Vert P_{p}^{\perp}\left[\nabla\cdot u\right]\right\Vert^{2} \lesssim \Vert u \Vert^{2}_{k+1} h^{2k}, \end{align} (6.11) and can therefore be cast as (6.7). A similar strategy can be followed for the corresponding terms in $$R_{11}(t)$$, $$R_{12}(t)$$ and $$R_{13}(t)$$. The bound for this last term requires some more work. The first term can be treated as above, using the step just described for the term with $${\mathscr{U}}$$ and the same strategy as in (6.10) for $$\hat{u}_{h}\cdot \nabla{\varSigma }$$. The second term in the bound of $$R_{13}(t)$$ can be absorbed by Q(t). It only remains to treat the last term, which is particularly important because it is one of the terms responsible for the need to have $$k \geqslant d/2$$. It can be bounded as follows: \begin{align} \left\Vert P_{\sigma}^{\perp}\left[g^{*}\left(\hat{u}_{h},F\right)\right]\right\Vert^{2} &\lesssim \left\Vert P_{\sigma}^{\perp}\left[g^{*}\left(\hat{u}_{h}-u,F\right)\right]\right\Vert^{2}+\left\Vert P_{\sigma}^{\perp}\left[g^{*}\left(u,F\right)\right]\right\Vert^{2}\nonumber\\ &\lesssim \left\Vert \nabla\left(\hat{u}_{h}-u\right)\right\Vert{}^{2}\left\Vert F\right\Vert{}_{\infty}^{2}+\left\Vert \nabla u\right\Vert{}_{\infty}^{2}\left\Vert F\right\Vert{}^{2}\nonumber\\ &\lesssim \left(h^{-d}\left\Vert \nabla\left(\hat{u}_{h}-u\right)\right\Vert{}^{2}+\left\Vert \nabla u\right\Vert{}_{\infty}^{2}\right)\left\Vert F\right\Vert{}^{2}, \end{align} (6.12) where an inverse inequality has been used in the last step. Since $$\hat{U}_{h}\in \mathscr{B}_{h}$$, $$\left \Vert \nabla \left (\hat{u}_{h}-u\right )\right \Vert^{2}$$ is of order $$h^{2k}$$, and the condition $$k\geqslant d/2$$ allows us to guarantee that the term in parenthesis remains bounded as h → 0. Therefore, (6.12) (multiplied by $$\varepsilon$$) can be absorbed by Q(t). This concludes the analysis of $$R_{13}(t)$$. In the bound for $$R_{14}(t)$$, the second and third terms in the RHS have already appeared when dealing with $$R_{13}(t)$$, and thus only the first term needs to be bounded. This can be done as follows: \begin{align} \left\Vert P_{\sigma}^{\perp}\left[g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right]\right\Vert^{2} &\lesssim \left\Vert g\left(\hat{u}_{h}-u,\hat{\sigma}_{h}-\sigma\right)\right\Vert{}^{2} +\left\Vert g\left(u,\hat{\sigma}_{h}-\sigma\right)\right\Vert{}^{2} \nonumber\\ &\quad +\left\Vert P_{\sigma}^{\perp}\left[ g\left(u,\sigma\right) \right] \right\Vert^{2}+\left\Vert g\left(\hat{u}_{h}-u,\sigma\right)\right\Vert{}^{2}, \end{align} (6.13) where the four terms of the RHS can be bounded using (2.4) and inverse inequalities. Using the fact that $$\hat{U}_{h}\in \mathscr{B}_{h}$$, $$k \geqslant d/2$$ and Assumption 6.1, it follows that all terms contributing to bound $$R_{14}(t)$$ (with the adequate factor) can be written as indicated in (6.6). Term $$R_{15}(t)$$ in (6.8) corresponds to the traction or rotational term of the Oldroyd-B constitutive model, and can be written as: \begin{align} R_{15}(t) & =-\frac{\lambda}{2\mu}(g\left(u,\sigma\right),F)+\frac{\lambda}{2\mu}(g\left(\hat{u}_{h},\hat{\sigma}_{h}\right),F)\nonumber\\ & = -\frac{\lambda}{2\mu}(g\left(u,\sigma-\hat{\sigma}_{h}\right),F)+\frac{\lambda}{2\mu}(g\left(u-\hat{u}_{h},\sigma-\hat{\sigma}_{h}\right),F) -\frac{\lambda}{2\mu}(g\left(u-\hat{u}_{h},\sigma\right),F). \end{align} (6.14) The three terms in the RHS can be bounded using Cauchy–Schwarz’s and Young’s inequalities and the inverse inequality (2.2). When these inequalities are used in (6.14) it is seen that we again recover (6.6). The remaining $$R_{i}$$ terms are all multiplied by $$\gamma$$. To distinguish these terms from the others we use $$\xi$$ instead of $$\varepsilon$$. The first of them can be bounded using Cauchy–Schwarz’s and Young’s inequalities, as well as the inverse inequality, as $$R_{16}(t) \leqslant (1-\beta)\gamma\alpha_{u}\left\Vert \nabla\varPi\right\Vert \left\Vert{P_{u,0}\left[\nabla G\right]}\right\Vert \lesssim \frac{1}{\mu}(1-\beta)\gamma\frac{1}{\xi}\left\Vert \varPi\right\Vert{}^{2} + (1-\beta)\gamma\mu h^{-2} {\alpha_{u}^{2}} \xi\left\Vert{P_{u,0}\left[\nabla G\right]}\right\Vert{}^{2}.$$ The next term requires some more elaboration. Observe first that $$( \partial_{t}(u-u_{h}),{P_{u,0}\left[\nabla G\right]}) = ( \partial_{t} u,{P_{u,0}\left[\nabla G\right]}) + ( \partial_{t} \nabla\cdot u_{h}, G).$$ Since u is divergence free and vanishes on $$\partial \varOmega$$, we have that $$( \partial_{t} u,{P_{u,0}\left[\nabla G\right]}) = ( \partial_{t} u, P_{u}[\nabla G] - \nabla G) = -\left(P_{u}^{\bot}[\partial_{t} u], P_{u}^{\bot}[\nabla G]\right).$$ Making use of the continuity equation (5.2), we can write $$( \partial_{t} \nabla\cdot u_{h}, G) = -\,\alpha_{u} \left(\partial_{t} P_{u}^{\bot} [\nabla p_{h}], P_{u}^{\bot} [\nabla G]\right) = \alpha_{u} \left(\partial_{t} P_{u}^{\bot} [\nabla G], P_{u}^{\bot} [\nabla G]\right) - \alpha_{u} \left(\partial_{t} P_{u}^{\bot} [\nabla{\mathscr{P}} ], P_{u}^{\bot} [\nabla G]\right).$$ From these two equalities, we obtain \begin{align*} R_{17}(t) =& \,(1-\beta)\gamma\alpha_{u}\rho ( \partial_{t}(u-u_{h}),{P_{u,0}\left[\nabla G\right]}) \\ =& -(1-\beta)\gamma\alpha_{u}\rho \left( P_{u}^{\bot}[\partial_{t} u], P_{u}^{\bot} [\nabla G]\right)\\ & + (1-\beta){\gamma\alpha_{u}^{2}} \rho \frac{1}{2}\frac{\mathrm{d}}{\mathrm{d} t} \left\Vert P_{u}^{\bot} [\nabla G] \right\Vert^{2} - (1-\beta){\gamma\alpha_{u}^{2}} \rho \left(\partial_{t} P_{u}^{\bot} [\nabla{\mathscr{P}}], P_{u}^{\bot} [\nabla G]\right). \end{align*} After using Schwarz’s and Young’s inequalities, the first term in the RHS can be treated as $$R_{4}(t)$$ (note that $$\varLambda = P_{u}^{\bot } [u]$$), and the third term can be treated as $$R_{11}(t)$$; therefore, they both can be cast as (6.7). It only remains to deal with the second term. Let $$\chi_{h} (G(t)) := (1-\beta)\alpha_{u} \left\Vert P_{u}^{\bot} [\nabla G(t)] \right\Vert^{2} \geqslant 0.$$ After integrating in time, in (6.9) we have the terms \begin{align} {\int_{0}^{T}} \chi_{h} (G(t))\,\mathrm{d} t + \dots \lesssim \alpha_{u} \rho \gamma \chi_{h}(G(T)) + \dots, \end{align} (6.15) the one on the RHS coming from $$R_{17}(t)$$. Since $$\hat{U}_{h}\in \mathscr{B}_{h}$$ and because of the boundedness of U described in Assumption 6.1, the RHS in (5.14) is bounded for all t. Since, as a consequence of the proof of Theorem 5.2, $$p_{h}$$ is continuous in $$[0,\infty )$$, we have that $$\alpha _{u} \Vert P_{u}^{\bot } [\nabla p_{h}(t)]\Vert ^{2}$$ is continuous and bounded in t. Assumption 6.1 implies that $$\alpha _{u} \Vert P_{u}^{\bot } [\nabla p(t)]\Vert ^{2}$$ is also continuous and bounded in t, and therefore $$\chi _{h} (G(T))$$ is continuous and bounded for all T. In this case, there exists $$T_{0}$$ such that \begin{align} \alpha_{u} \rho \chi_{h}(G(T)) \leqslant{\int_{0}^{T}} \chi_{h}(G(t))\,\mathrm{d}t =:\mathscr{G}(T), \quad \textrm{for all}~T \geqslant T_{0}. \end{align} (6.16) Let us prove this. The LHS is bounded in T, and thus the result is obvious if $$\mathscr{G}(T)$$ is unbounded. Suppose now that $$\mathscr{G}(T) \nearrow M < \infty$$ as $$T \to \infty$$. In this case, $$\chi _{h}(G(t)) \to 0^{+}$$ as $$t\ \to\ \infty$$. Thefore, there must exist $$T_{1}$$ such that $$\mathscr{G}(T)> M/2$$ for all $$T> T_{1}$$ and $$T_{0}> T_{1}$$ such that $$\alpha _{u} \rho \chi _{h}(G(T)) < M/2$$ for all $$T> T_{0}$$, so that (6.16) holds. Once (6.16) is proved, it is observed that the term in the RHS of (6.15) can be absorbed by the one in the LHS for an appropriate $$\gamma$$. This concludes the analysis of $$R_{17}(t)$$. The following terms can be controlled using the inverse inequality and Young’s inequality: \begin{align*} R_{18}(t) & \lesssim \beta(1-\beta)\gamma\alpha_{u}\mu\left(\left\Vert \nabla^{s}E\right\Vert +\left\Vert \nabla^{s}\varLambda\right\Vert \right)\left\Vert \nabla^{s}\left({P_{u,0}\left[\nabla G\right]}\right)\right\Vert\!,\\ R_{19}(t) & \leqslant (1-\beta)\gamma\alpha_{u}\left(\left\Vert F\right\Vert +\left\Vert \varGamma\right\Vert \right)\left\Vert \nabla^{s}\left({P_{u,0}\left[\nabla G\right]}\right)\right\Vert\!. \end{align*} To see that the bounds obtained fall within the structure (6.6), we have to use the expression of the stabilization parameter $$\alpha _{u}$$ to check that $$\mu h^{-2}\lesssim \alpha _{u}^{-1}$$ and take $$\gamma$$ small enough to balance $$\xi ^{-1}$$, which may need to be large. The next terms to consider come from the stabilization \begin{align*} & R_{20}(t) =\!(1-\beta){\gamma\alpha_{u}^{2}}\left(P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla\left(u_{h}-{\mathscr{U}} +{\mathscr{U}}\right)\right],P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla\left({P_{u,0}\left[\nabla G\right]}\right)\right]\right)\!,\\ & R_{21}(t) =\!(1-\beta)\gamma\alpha_{u}\alpha_{p}\left(P_{p}^{\perp}\left[\nabla\cdot\left(u_{h}-{\mathscr{U}} +{\mathscr{U}}\right)\right],P_{p}^{\perp}\left[\nabla\cdot\left({P_{u,0}\left[\nabla G\right]}\right)\right]\right)\!,\\ & R_{22}(t) =\!(1\!-\!\beta)\gamma\alpha_{u}\alpha_{\sigma}\!\left(\!P_{\sigma}^{\perp}\left[\!(1\!-\!\beta)\nabla^{s}\!\left(u_{h}\!-{\mathscr{U}} \!+{\mathscr{U}}\right)\!-\!\frac{\lambda}{2\mu}\hat{u}_{h}\cdot\nabla\left(\sigma_{h}\!-\!{\varSigma} \!+\!{\varSigma}\right)\!\right]\!, P_{\sigma}^{\perp}\left[\nabla^{s}\!\left({P_{u,0}\left[\nabla G\right]}\right)\right]\right)\!,\\ & R_{23}(t) =\!(1-\beta)\gamma\alpha_{u}\alpha_{\sigma}\left(P_{\sigma}^{\perp}\left[\frac{\lambda}{2\mu}g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right],P_{\sigma}^{\perp}\left[\nabla^{s}\left({P_{u,0}\left[\nabla G\right]}\right)\right]\right). \end{align*} All these terms can be bounded using a similar strategy as before. With this we conclude the analysis of the last term in (6.8) (or in (6.9)). The second term, S(t), is easy to treat; after using Young’s inequality, the first term that arises can be absorbed by Q(t) and the second bounded as in (6.12). We have written S(t) independently to emphasize that it arises because of the linearization of the rotational term: since we evaluate it with $$\hat{u}_{h}$$ and $$\hat{\sigma }_{h}$$, it does not contribute to stability (that is to say, its contribution does not appear in Q(t)). We could however have considered it as one more $$R_{j}(t)$$ term. It only remains to deal with the convective terms N(t) in (6.8). The first two of these terms can be written as follows: \begin{align} & \rho(1-\beta)\left(\tilde{c}\left(u,u,E\right)-\tilde{c}\left(\hat{u}_{h},u_{h},E\right)\right)\nonumber \\ & \quad = \rho(1-\beta)\left(\tilde{c}\left(\left(u-\hat{u}_{h}\right),u,E\right)-\tilde{c}\left(\left(u-\hat{u}_{h}\right),\varLambda,E\right)+\tilde{c}\left(u,\varLambda,E\right)\right)\!. \end{align} (6.17) The three RHS terms of the above equation can be easily bounded fitting the structure (6.6)–(6.7). A similar procedure can be applied to the convective terms arising from the constitutive equation, which can be written as: $$\frac{\lambda}{2\mu}\left(\tilde{c}\left(u,\sigma,F\right)-\tilde{c}\left(\hat{u}_{h},\sigma_{h},F\right)\right) = \frac{\lambda}{2\mu}\left(\tilde{c}\left(\left(u-\hat{u}_{h}\right),\sigma,F\right)-\tilde{c}\left(\left(u-\hat{u}_{h}\right),\varGamma,F\right)+\tilde{c}\left(u,\varGamma,F\right)\right)\!.$$ The remaining convective terms to be controlled in the LHS of (6.8) can be written as \begin{align*} &(1-\beta)\gamma\alpha_{u}\rho\left(\tilde{c}\left(u,u,{P_{u,0}\left[\nabla G\right]}\right)-\tilde{c}\left(\hat{u}_{h},u_{h},{P_{u,0}\left[\nabla G\right]}\right)\right)\\ & \quad = (1-\beta)\gamma\alpha_{u}\rho\left( \tilde{c}(u-\hat{u}_{h},u,{P_{u,0}\left[\nabla G\right]}) - \tilde{c}(u-\hat{u}_{h},u -u_{h},{P_{u,0}\left[\nabla G\right]}) + \tilde{c}(u,u-u_{h},{P_{u,0}\left[\nabla G\right]}) \right)\!. \end{align*} The strategy to bound these terms is similar to what we have used heretofore. This concludes the proof of the third step. Step 4. Fixed point theorem. According to Step 3, the ball $$\mathscr{B}_{h}$$ is invariant under the map $$\delta$$, that is to say, $$\delta \left (\mathscr{B}_{h}\right )\subset \mathscr{B}_{h}$$. Therefore, applying Schauder’s fixed point theorem, we conclude that there exists $$(u_{h},p_{h},\sigma _{h})\in \mathscr{B}_{h}$$ such that $$(u_{h},p_{h},\sigma _{h}) = \delta (u_{h},p_{h},\sigma _{h})$$, which, in view of the definition of $$\delta$$, implies that $$(u_{h},p_{h},\sigma _{h})$$ is a solution of (4.1). Because of the definition of $$\mathscr{B}_{h}$$ in (6.2), we also obtain an error estimate. 7. Conclusions In this article we have presented the numerical analysis of a stabilized finite element approximation proposed to solve viscoelastic fluid flows, in the nonlinear time-dependent case. This analysis has confirmed the numerical results obtained in other works, where the method was proposed and tested in nonlinear examples. In particular, we have shown this using a fixed point theorem under suitable regularity conditions in velocity, stress and pressure. As it is usual for nonlinear problems involving the Navier–Stokes equations, the estimates do not show information about their behaviour when the local Reynolds and Weissenberg numbers increase. On the other hand, some control over the pressure, although very weak, has been obtained. Funding E.C. thanks the support given by the Chilean Council for Scientific and Technological Research (CONICYT-FONDECYT 11160160), and from the Universidad de S. de Chile. R.C. gratefully acknowledges the support received from the Catalan Government through the Catalan Institution for Research and Advanced Studies (ICREA) Acadèmia Research Program. References Baaijens , F . ( 1998 ) Mixed finite element methods for viscoelastic flow analysis: a review . J. Nonnewton. Fluid Mech. , 79 , 361 – 385 . Google Scholar CrossRef Search ADS Badia , S. & Codina , R. ( 2009 ) On a multiscale approach to the transient Stokes problem: dynamic subscales and anisotropic space–time discretization . Appl. Math. Comput. , 207 , 415 – 433 . Badia , S. , Codina , R. & Gutiérrez-Santacreu , J. V. ( 2010 ) Long-term stability estimates and existence of a global attractor in a finite element approximation of the Navier–Stokes equations with numerical subgrid scale modeling . SIAM J. Numer. Anal. , 48 , 1013 – 1037 . Google Scholar CrossRef Search ADS Baranger , J. & Sandri , D. ( 1992 ) Finite element approximation of viscoelastic fluid flow: existence of approximate solutions and error bounds . Numer. Math. , 63 , 13 – 27 . Google Scholar CrossRef Search ADS Baranger , J. & Wardi , S. ( 1995 ) Numerical analysis of a FEM for a transient viscoelastic flow . Comput. Methods Appl. Mech. Eng. , 125 , 171 – 185 . Google Scholar CrossRef Search ADS Barrett , J. & Boyaval , S. ( 2011 ) Existence and approximation of a (regularized) oldroyd-b model . Math. Models Meth. Appl. Sci. , 21 , 1783 – 1837 . Google Scholar CrossRef Search ADS Bird , R. B. , Armstrong , R. C. & Hassager , O. ( 1987a ) Dynamics of Polymeric Liquids, vol. 1: Fluid Mechanics , 2nd edn. New York : Wiley . Bird , R. B. , Armstrong , R. C. & Hassager , O. ( 1987b ) Dynamics of Polymeric Liquids, vol. 2: Kinetic Theory , 2nd edn. New York : Wiley . Bird , R. B. , Stewart , W. E. & Lightfoot , E. N. ( 2002 ) Transport Phenomena . New York : Wiley . Bogaerds , A. , Verbeeten , W. , Peters , G. & Baaijens , F. ( 1999 ) 3D viscoelastic analysis of a polymer solution in a complex flow . Comput. Meth. Appl. Mech. Eng. , 180 , 413 – 430 . Google Scholar CrossRef Search ADS Bonito , A. & Burman , E. ( 2007 ) A continuous interior penalty method for viscoelastic flows . SIAM J. Sci. Comput. , 30 , 1156 – 1177 . Google Scholar CrossRef Search ADS Bonito , A. , Clément , P. & Picasso , M. ( 2007 ) Mathematical and numerical analysis of a simplified time-dependent viscoelastic flow . Num. Math. , 107 , 213 – 255 . Google Scholar CrossRef Search ADS Burman , E. & Fernández , M. ( 2007 ) Continuous interior penalty finite element method for the time-dependent Navier–Stokes equations: space discretization and convergence . Num. Math. , 107 , 39 – 77 . Google Scholar CrossRef Search ADS Castillo , E. , Baiges , J. & Codina , R. ( 2015 ) Approximation of the two-fluid flow problem for viscoelastic fluids using the level set method and pressure enriched finite element shape functions . J. Nonnewton. Fluid Mech. , 225 , 37 – 53 . Google Scholar CrossRef Search ADS Castillo , E. & Codina , R. ( 2014a ) Stabilized stress-velocity-pressure finite element formulations of the Navier-Stokes problem for fluids with non-linear viscosity . Comput. Methods Appl. Mech. Eng. , 279 , 554 – 578 . Google Scholar CrossRef Search ADS Castillo , E. & Codina , R. ( 2014b ) Variational multi-scale stabilized formulations for the stationary three-field incompressible viscoelastic flow problem . Comput. Methods Appl. Mech. Eng. , 279 , 579 – 605 . Google Scholar CrossRef Search ADS Castillo , E. & Codina , R. ( 2015 ) First, second and third order fractional step methods for the three-field viscoelastic flow problem . J. Comput. Phys. , 296 , 113 – 137 . Google Scholar CrossRef Search ADS Castillo , E. & Codina , R. ( 2017a ) Finite element approximation of the viscoelastic flow problem: a non-residual based stabilized formulation . Comput. Fluids , 142 , 72 – 78 . Google Scholar CrossRef Search ADS Castillo , E. & Codina , R. ( 2017b ) Numerical analysis of a stabilized finite element approximation for the three-field linearized viscoelastic fluid problem using arbitrary interpolations . ESAIM Math. Model. Num. Anal. , 51 , 1407 – 1427 . Ciarlet , P. ( 2013 ) Linear and Nonlinear Functional Analysis with Applications . Philadelphia, PA, USA : Society for Industrial and Applied Mathematics . Codina , R . ( 2000 ) Stabilization of incompressibility and convection through orthogonal sub-scales in finite element methods . Comput. Methods Appl. Mech. Eng. , 190 , 1579 – 1599 . Google Scholar CrossRef Search ADS Codina , R . ( 2009 ) Finite element approximation of the three-field formulation of the Stokes problem using arbitrary interpolations . SIAM J. Numer. Anal. , 47 , 699 – 718 . Google Scholar CrossRef Search ADS Codina , R. & Blasco , J. ( 2000 ) Stabilized finite element method for the transient Navier–Stokes equations based on a pressure gradient projection . Comput. Methods Appl. Mech. Eng. , 182 , 277 – 300 . Google Scholar CrossRef Search ADS Codina , R. , Principe , J. , Guasch , O. & Badia , S. ( 2007 ) Time dependent subscales in the stabilized finite element approximation of incompressible flow problems . Comput. Methods Appl. Mech. Eng. , 196 , 2413 – 2430 . Google Scholar CrossRef Search ADS Ern , A. & Guermond , J.-L. ( 2004 ) Theory and Practice of Finite Elements , vol. 159, 1st edn. New York : Springer . Google Scholar CrossRef Search ADS Ervin , V. , Lee , H. & Ntasin , L. ( 2005 ) Analysis of the Oseen-viscoelastic fluid flow problem . J. Nonnewton. Fluid Mech. , 127 , 157 – 168 . Google Scholar CrossRef Search ADS Ervin , V. & Miles , W. ( 2003 ) Approximation of time-dependent viscoelastic fluid flow: SUPG approximation . SIAM J. Numer. Anal. , 41 , 457 – 486 . Google Scholar CrossRef Search ADS Ervin , V. & Miles , W. ( 2005 ) Approximation of time-dependent, multi-component, viscoelastic fluid flow . Comput. Methods Appl. Mech. Eng. , 194 , 2229 – 2255 . Google Scholar CrossRef Search ADS Fernández-Cara , E. , Guillén , F. & Ortega , R. ( 2002 ) Mathematical modeling and analysis of viscoelastic fluids of the Oldroyd kind . Handbook of Numerical Analysis , vol. 8. Amsterdam : North-Holland , pp. 543 – 661 . Fortin , M. & Fortin , A. ( 1989 ) A new approach for the FEM simulation of viscoelastic flows . J. Nonnewton. Fluid Mech. , 32 , 295 – 310 . Google Scholar CrossRef Search ADS Girault , V. & Raviart , P.-A. ( 1986 ) Finite Elements for the Navier–Stokes Equations . Berlin : Springer . Google Scholar CrossRef Search ADS Guillopé , C. & Saut , J. ( 1990 ) Existence results for the flow of viscoelastic fluids with a differential constitutive law . Nonlinear Anal. Theory Methods Appl. , 15 , 849 – 869 . Google Scholar CrossRef Search ADS Howell , J . ( 2009 ) Dual-mixed finite element approximation of Stokes and nonlinear Stokes problems using trace-free velocity gradients . J. Comput. Appl. Math. , 231 , 780 – 792 . Google Scholar CrossRef Search ADS Hughes , T. , Feijóo , G. , Mazzei , L. & Quincy , J.-B. ( 1998 ) The variational multiscale method-a paradigm for computational mechanics . Comput. Methods Appl. Mech. Eng. , 166 , 3 – 24 . Google Scholar CrossRef Search ADS Lions , P. L. & Masmoudi , N. ( 2000 ) Global solutions for some Oldroyd models of non-Newtonian flows . Chin. Ann. Math. Ser. B , 21 , 131 – 146 . Google Scholar CrossRef Search ADS Marchal , J. & Crochet , M. ( 1987 ) A new mixed finite element for calculating viscoelastic flow . J. Nonnewton. Fluid Mech. , 26 , 77 – 114 . Google Scholar CrossRef Search ADS Pani , A. K. & Yuan , J. Y. ( 2005 ) Semidiscrete finite element Galerkin approximations to the equations of motion arising in the Oldroyd model . IMA J. Numer. Anal. , 25 , 750 – 782 . Google Scholar CrossRef Search ADS Picasso , M. & Rappaz , J. ( 2001 ) Existence, a priori and a posteriori error estimates for a nonlinear three-field problem arising from Oldroyd-B viscoelastic flows . ESAIM: Math. Model. Numer. Anal. , 35 , 879 – 897 . Google Scholar CrossRef Search ADS Renardy , M . ( 1985 ) Existence of slow steady flows of viscoelastic fluids with differential constitutive equations . ZAMM/J. Appl. Math. Mech. , 65 , 449 – 451 . Google Scholar CrossRef Search ADS Renardy , M. ( 2000 ) Mathematical analysis of viscoelastic flows . CBMS-NSF Regional Conference Series in Applied Mathematics . SIAM : Philadelphia . Google Scholar CrossRef Search ADS Sandri , D . ( 1994 ) Finite element approximation of viscoelastic fluid flow: existence of approximate solutions and error bounds. Continuous approximation of the stress . SIAM J. Numer. Anal. , 31 , 362 – 377 . Google Scholar CrossRef Search ADS Sureshkumar , R. & Beris , A. ( 1995 ) Effect of artificial stress diffusivity on the stability of numerical calculations and the flow dynamics of time-dependent viscoelastic flows . J. Nonnewton. Fluid Mech. , 60 , 53 – 80 . Google Scholar CrossRef Search ADS © The Author(s) 2018. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) For permissions, please e-mail: journals. permissions@oup.com http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

# Time-dependent semidiscrete analysis of the viscoelastic fluid flow problem using a variational multiscale stabilized formulation

, Volume Advance Article – Apr 24, 2018
28 pages

/lp/ou_press/time-dependent-semidiscrete-analysis-of-the-viscoelastic-fluid-flow-0xLi2dGifV
Publisher
Oxford University Press
ISSN
0272-4979
eISSN
1464-3642
D.O.I.
10.1093/imanum/dry018
Publisher site
See Article on Publisher Site

### Abstract

Abstract In this article we analyse a stabilized finite element formulation recently proposed to approximate viscoelastic fluid flows. The formulation has shown to have accuracy and robustness in the different benchmarks tested in the viscoelastic framework, and permitting the use of equal interpolation of the unknown fields. We first present results about a linearized subproblem, for which well-posedness and stability results can be proved. Then, the semidiscrete nonlinear time-dependent case is addressed using a fixed point theorem, which allows us to prove existence of a semidiscrete solution, along with error estimates. 1. Introduction Viscoelastic fluids are a specific type of nonNewtonian fluids. They are characterized by having complex and high molecular weight molecules with many internal degrees of freedom (Bird et al., 2002). The classical examples of this type of fluids are the polymer solutions and molten polymers. The basic feature of polymeric fluids is the presence of long-chain molecules. In a flow, these chain molecules are stretched out by the drag forces exerted on them by the surrounding fluid (Renardy, 2000). The natural tendency of the molecule to retract from this stretched configuration generates an elastic force, which contributes to the macroscopic stress tensor, and for this reason they are called viscoelastic fluids. The interest for fluids of this kind has increased in recent years, due to the connections with the industrial applications. This motivates the numerical and mathematical analysis of the governing equations (see, e.g., Fernández-Cara et al., 2002). For viscoelastic fluid flows, in contrast to the Navier–Stokes equations, well-posedness for general models is not well understood. For initial value problems, the existence of solutions has been proved only locally in time. Global existence in time of solutions has been proved only if the initial conditions are small perturbations of the rest state, and for the steady state case existence of solutions can be proved only for small perturbations of the Newtonian case (see Renardy, 2000 for a comprehensive review). The existence of slow steady flows of viscoelastic fluids using differential constitutive equations was proved in the study by Renardy (1985) for Hilbert spaces. In this work, the author used an iterative method to show that the solution can be bounded in a certain norm, and then he proved that all iterates converge in a weaker norm. For the time-dependent case, existence of solutions locally in time, and for small data globally in time, has been proved for Hilbert spaces in the study by Guillopé & Saut (1990). The extension to Banach spaces and a complete review of uniqueness, regularity, well-posedness and stability results can be found in Fernández-Cara et al. (2002). The existence of global weak solutions for general initial conditions using a corotational Oldroyd-B model has been proved in the study by Lions & Masmoudi (2000) using a simplification (without physical justification), which consists in replacing the velocity gradient in the stress equation by its skew-symmetric part. In the study by Barrett & Boyaval (2011) the authors proved global existence of weak solutions in two dimensions to the Oldroyd-B model regularized with the introduction of a diffusion term in the constitutive equation, and assuming homogeneous natural boundary conditions associated to this term. The proof in this paper is based on a mixed finite element interpolation of the problem. The introduction of the stress diffusion can be physically justified. An analysis of the effects it has on the numerical approximation can be found in the study by Sureshkumar & Beris (1995). From a finite element perspective, the finite element approximation of the flow of viscoelastic fluids presents several numerical difficulties. On the one hand, there are all the problems inherited from the incompressible Navier–Stokes equations, mainly the compatibility between the velocity–pressure approximation and the treatment of the nonlinear advective term. But, on top of that, now the constitutive equation is highly nonlinear, with an advective term that may lead to both global and local oscillations in the numerical approximation. Moreover, even in the case of smooth solutions, it is necessary to meet some additional compatibility conditions between the velocity and the stress interpolation in order to control velocity gradients. Elements that satisfy the compatibility requirements velocity–pressure and stress–velocity are scarce, particularly in the three-dimensional case (see, e.g., Marchal & Crochet, 1987; Fortin & Fortin, 1989; Bogaerds et al., 1999). In the study by Baaijens (1998) one can find a good review of mixed methods that satisfy the two required compatibility conditions. In the finite element framework, the work of Baranger & Sandri (1992) was one of the first where the existence of an approximate solution and error bounds were given for an Oldroyd-B fluid, using Brouwer’s fixed-point theorem, discontinuous interpolation for the elastic stresses and the method of Lesaint–Raviart for the convection of the extra stress tensor in the stationary case. Later, Sandri extended (Sandri, 1994) the study to continuous approximation of the stress field, using $$P_{1}$$-$$P_{2}$$-$$P_{1}$$ interpolation (linear-quadratic-linear) for $$\sigma$$ (stress), u (velocity) and p (pressure), respectively, and the Streamline-Upwind/Petrov-Galerkin (SUPG) method to treat the convective term in the constitutive equation. The time-dependent case of the same continuous interpolation was analysed in the study by Baranger & Wardi (1995). In the same finite element context, Picasso & Rappaz (2001) analysed a stationary nonconvective Oldroyd-B problem proving a priori and a posteriori error estimates. In this work, the authors used the Galerkin Least Squares method to stabilize the momentum equation and the Elastic Viscous Split Stress (EVSS) scheme for the constitutive equation. The extension to the time-dependent case was treated in the study by Bonito et al. (2007) for the same simplified Oldroyd-B problem, proving global existence in time in Banach spaces for small data. For a Stokes/Oldroyd-B linearized problem, Bonito & Burman (2007) presented optimal a priori error estimates using the Interior Penalty method. A similar problem was studied in the study by Ervin et al. (2005) for the steady state case, but using the Johnson–Segalman linearized constitutive model, proving existence and uniqueness of the continuous solution and of a finite element approximation under a small data assumption. In another work, Ervin & Miles (2003) analysed the Oldroyd-B time-dependent case, both in the semidiscrete and in the fully discrete cases using the SUPG method, and proving existence of a solution and deriving a priori error estimates for the numerical approximation, assuming a Taylor–Hood pair approximation for the velocity and pressure and a continuous approximation for the viscoelastic stresses. The same authors extended later the analysis to a two-fluid flow problem in the study by Ervin & Miles (2005), giving a priori error estimates for the approximation in terms of the mesh and time discretiation parameters. In the study by Pani & Yuan (2005) the authors analysed the time behaviour of the viscoelastic Oldroyd model in two dimensions using a Galerkin formulation in space; in this work, the stress is eliminated through a proper projection operator, resulting in an integro-differential equation in terms of velocity and pressure. The stabilized finite element formulation analysed in this work has its roots in the Variational Multiscale (VMS) method introduced in the study by Hughes et al. (1998) for the scalar convection–diffusion-reaction problem, and extended later to the Stokes problem in the study by Codina (2000), where the space of the subgrid scales is taken orthogonal to the finite element space. As we shall see, this is an important ingredient in the method analysed, which consists in a sort of orthogonal term-by-term stabilized formulation. The key idea behind a VMS method consists in splitting the unknowns of the problem in two scales, the finite element one and the unresolvable one, called subgrid scale. The latter needs to be approximated in a simple manner in terms of the former, so as to capture its main effect and yield a globally stable formulation for the finite element unknown, keeping therefore the number of degrees of freedom of the Galerkin method. The objective of this paper is to analyse numerically the stabilized finite element formulation proposed in the study by Castillo & Codina (2014b) for the time-dependent viscoelastic flow problem. This formulation has shown to have very good accuracy and robustness in stationary (Castillo & Codina, 2014b) and time-dependent (Castillo & Codina, 2015; Castillo et al., 2015) cases. The numerical analysis of a linearized stationary case was performed in the study by Castillo & Codina (2017b), where in addition jump functions were added to permit arbitrary discontinuous interpolations of pressure and stresses. As it is usual in the analysis of numerical approximations to flow problems, even in the Newtonian case, our analysis is based on some stringent regularity assumptions on the continuous solution. Apart from analysing a nonstandard stabilized formulation, the novelty of this work is also the treatment of some of the terms that appear in the analysis. The work is organized as follows. Section 2 defines some notation and presents general results used in the subsequent analysis. In Section 3 we present the problem to be solved and its Galerkin finite element approximation, explaining the sources of the numerical instability. Section 4 contains a description of the stabilized finite element formulation analysed. Section 5 is devoted to the numerical analysis of a linearized time-dependent subproblem, where the stability of the method and the existence and uniqueness of the solution in the semidiscrete linearized case are proved. Section 6 analyses the nonlinear case, where existence of a solution is proved using a fixed point theorem. Finally, conclusions and some remarks are summarized in Section 7. 2. Notation and preliminaries 2.1 Notation Let us introduce some notation used hereafter. As usual, given a domain $$\omega$$ of $$\mathbb{R}^{d}$$, d = 2, 3, $$W^{m,p}(\omega )$$ denotes the Sobolev space of functions whose distributional derivatives of order up to $$m\in \mathbb{N}$$ belong to $$L^{p}(\omega )$$, $$p\geqslant 1$$, endowed with the standard norm. For p = 2 we write $$W^{m,2}(\omega )\equiv H^{m}(\omega )$$. The space $${H_{0}^{1}}\left (\omega \right )$$ consists of functions in $$H^{1}\left ({\omega} \right )$$ vanishing on $$\partial \omega$$. The topological dual of $${H_{0}^{1}}\left (\omega \right )$$ is denoted by $$H^{-1}\left (\omega \right )$$, the corresponding duality pairing by $$\left \langle \cdot ,\cdot \right \rangle _{\omega }$$, and the $$L^{2}$$ inner product in $$\omega$$ (for scalars, vectors and tensors) is denoted by $$\left (\cdot ,\cdot \right )_{\omega }$$. The symbol $$\left \langle \cdot ,\cdot \right \rangle _{\omega }$$ is also used for the integral over $$\omega$$ of the product of two functions, whenever it makes sense. When $$\omega =\varOmega$$, the domain where the problem is posed, the subscript or the domain information are omitted. Referring to the norms used in the subsequent analysis, $$\left \Vert \cdot \right \Vert$$ represents the $$L^{2}\left (\varOmega \right )$$-norm, $$\left \Vert \cdot \right \Vert_{L^{p}}$$ the $$L^{p}\left (\varOmega \right )$$-norm ($$2< p <\infty$$), $$\left \Vert \cdot \right \Vert_{\infty }$$ the $$L^{\infty }\left (\varOmega \right )$$-norm, $$\left \Vert \cdot \right \Vert_{k}$$ the $$H^{k}\left (\varOmega \right )$$-norm, $$\left \Vert \cdot \right \Vert_{l,p}$$ the norm for the space $$W^{l,p}\left (\varOmega \right )$$ (in particular, $$\left \Vert \cdot \right \Vert_{1,\infty }$$ is the $$W^{1,\infty }\left (\varOmega \right )$$-norm) and $$\left \vert \cdot \right \vert_{l,p}$$ the seminorm. The symbol $$\lesssim$$ will be used for $$\leqslant$$ up to constants independent of the physical parameters and the parameters of the numerical discretization. Given a vector field v, its symmetrical gradient will be denoted by $$\nabla ^{s}v$$; it is defined as $$\nabla^{s}v=\frac{1}{2}\left(\nabla v+\left(\nabla v\right)^{T}\right)\!.$$ Finally, the temporal derivative will be written as $$\partial _{t}$$. 2.2 Preliminaries Let us introduce some basic results. The gradient of a vector field v can be bounded in terms of its symmetrical gradient as in Ciarlet (2013): \begin{align} \left\Vert \nabla v\right\Vert \leqslant\sqrt{2}\left\Vert \nabla^{s} v\right\Vert,\quad v\in\left({H_{0}^{1}}\left(\varOmega\right)\right)^{d}\!.\nonumber \end{align} Poincaré–Friedrich’s and Korn’s inequalities read as: for $$v\in \left ({H_{0}^{1}}\left (\varOmega \right )\right )^{d}$$, there exists constants $$c_{\mathrm{PF}}$$ and $$c_{\textrm{K}}$$ such that \begin{align} \left\Vert v\right\Vert \leqslant c_{\mathrm{PF}}\left\Vert \nabla v\right\Vert,\;\;\;\left\Vert v\right\Vert{}_{1} \leqslant c_{\mathrm{K}}\left\Vert \nabla^{s} v\right\Vert\!.\nonumber \end{align} From the Sobolev embedding theorems (see Girault & Raviart, 1986, for example), we can bound the $$L^{4}(\varOmega )$$-norm in terms of the $$H^{1}(\varOmega )$$-norm as follows: \begin{align} \left\Vert v\right\Vert{}_{L^{4}}\leqslant C_{04}^{12} \left\Vert v\right\Vert{}_{1}. \end{align} (2.1) Constants $$c_{\mathrm{PF}}$$, $$c_{\mathrm{K}}$$ and $$C_{04}^{12}$$ depend on the shape and size of the domain $$\varOmega$$. We will use $${{c}_{\mathrm{PF}}}$$ as the constant to bound the $$L^{2}(\varOmega )$$-norm of a function by its $$H^{1}(\varOmega )$$-norm. For the finite element formulation that we shall consider, the discrete velocity field will not be pointwise divergence free and the use of the skew-symmetric counterpart of the convective term simplifies the subsequent analysis. Then, for u, v, $${w}\in ({H^{1}_{0}}(\varOmega ))^{d}$$ we define $$\tilde{c}\left({w},u, v\right)= \frac{1}{2}\left(c\left({w},u, v\right)-c\left({w}, v,u\right)\right),\quad \textrm{with}\quad c\left({w},u, v\right)= \langle{w}\cdot\nabla u, v\rangle.$$ The following properties are satisfied by this skew-symmetric form: $$\tilde{c}\left ({w},u, v\right )=c\left ({w},u, v\right )$$ when ∇⋅ w = 0 and either w = 0 or u = 0 or v = 0 on $$\partial \varOmega$$. $$\tilde{c}({w},u,u)=0$$, $${w},u\in (H^{1}(\varOmega ) )^{d}$$. $$\tilde{c}({w},u, v) \leqslant \left (C_{04}^{12}\right )^{2} \Vert{w}\Vert _{1} \Vert u\Vert _{1}\Vert v\Vert _{1}$$, $${w},u, v\in (H^{1}(\varOmega ))^{d}$$. $$\tilde{c}\left ({w},u, v\right ) \leqslant \Vert{w}\Vert _{\infty } ( \Vert u\Vert _{1}\Vert v\Vert + \Vert u\Vert \Vert v\Vert _{1})$$, $${w}\in (L^{\infty }(\varOmega ))^{d}$$, $$u, v\in (H^{1}(\varOmega ))^{d}$$. $$\tilde{c}\left ({w},u, v\right ) \leqslant \Vert v\Vert ( \Vert u\Vert _{\infty }\Vert \nabla{w}\Vert + \Vert \nabla u\Vert _{\infty }\Vert{w}\Vert )$$, $${w}\in ({H_{0}^{1}}(\varOmega ))^{d}$$, $$u\in (W^{1,\infty }(\varOmega ))^{d}$$, $$v\in (L^{2}(\varOmega ))^{d}$$. $$\tilde{c}\left ({w},u, v\right ) \leqslant \Vert v\Vert ( \Vert u\Vert \Vert \nabla{w}\Vert _{\infty } + \Vert \nabla u\Vert \Vert{w}\Vert _{\infty })$$, $${w}\in (W^{1,\infty }(\varOmega ))^{d}$$, $$u\in ({H^{1}_{0}}(\varOmega ))^{d}$$, $$v\in (L^{2}(\varOmega ))^{d}$$. Integration by parts has to be used to prove the last two properties. The same definition of the trilinear forms c and $$\tilde{c}$$ and with the same properties can be introduced when the last two arguments are tensor valued and with the appropriate regularity. Let us consider a finite element partition $$\mathscr{T}_{h} = \{ K\}$$ of the computational domain $$\varOmega$$. The diameter of an element domain $$K\in \mathscr{T}_{h}$$ is denoted by $$h_{K}$$ and the diameter of the partition, or mesh size, is defined as $$h=\max \left \{ h_{K}\mid K\in \mathscr{T}_{h}\right \}$$. We will consider for simplicity quasi-uniform families of meshes, and thus all the element diameters can be bounded above and below by constants multiplying h. Defining $$V_{h}:=\left \{ v_{h} : \varOmega \rightarrow \mathbb{R}~\vert ~ v_{h}\vert _{K}\in P_{k}\left (K\right )\:\forall \ K\in \mathscr{T}_{h}\right \}$$, $$P_{k} (K)$$ being the set of polynomials of degree k in K, we can write the following inverse inequality: there is a constant $$c_{\mathrm{inv}}$$, independent of the mesh size h, such that \begin{align} \left\Vert v_{h}\right\Vert{}_{l,\,p,\,K} \leqslant c_{\mathrm{inv}}h^{m-l+\min\left(0,\frac{d}{p}-\frac{d}{q}\right)}\left\Vert v_{h}\right\Vert{}_{m,\,q,\,K}, \end{align} (2.2) for all $$l\geqslant m, l\in [1,+\infty ]$$, and all finite element functions $$v_{h}$$ defined on $$K\in \mathscr{T}_{h}$$ (see Ern & Guermond, 2004, for example), where $$\left \Vert v_{h}\right \Vert_{l,\,p,\,K}$$ is the norm of $$v_{h}$$ in $$W^{l,\,p}(K)$$. For $$v\in W^{l,\,p}\left (\varOmega \right )$$ we have the following interpolation estimate: there exists a constant C independent of h, such that for all $$0\leqslant l\leqslant k+1$$, $$1\leqslant p\leqslant \infty$$, there holds \begin{align} \left\Vert v-P\left[v\right]\right\Vert{}_{L\,^{p}}\leqslant Ch^{l}\left| v \right|{}_{l,\,p}, \end{align} (2.3) where $$P\left [ v \right ]\in V_{h}$$ is the $$L^{2}\left (\varOmega \right )$$-projection of v in $$V_{h}$$ (see Ern & Guermond, 2004 for more details). From Sobolev’s embedding and the stability of P on finite element spaces we also have that \begin{align} \left\Vert \nabla\left(v-P\left[v\right]\right)\right\Vert{}_{\infty} \leqslant C_{\infty,3}\left\Vert v \right\Vert_{3}\!, \end{align} (2.4) where $$C_{\infty ,3}$$ is a positive constant independent of h. 3. Problem statement and Galerkin finite element discretization 3.1 The boundary value problem Let $$\varOmega$$ be an open set of $$\mathbb{R}^{d}$$ occupied by the fluid, assumed to be bounded and polyhedral, and let $$\partial \varOmega$$ be its boundary. Additionally, consider the time interval $$\left ]0,T\right [$$, with $$T<\infty$$. The incompressible and isothermal viscoelastic fluid flow problem can be written as $$\rho\left({\partial_{t}u}+u\cdot\nabla u\right)-\nabla\cdot\left(2\beta\mu\nabla^{s}u+\sigma\right)+\nabla p =f\qquad\textrm{in}\:\varOmega,\:t\in\left]0,T\right[,$$ (3.1) \begin{align} \nabla\cdot u &=0\qquad\textrm{in}\:\varOmega,\:t\in\left]0,T\right[, \end{align} (3.2) \begin{align} \frac{1}{2\mu}\sigma-(1-\beta)\nabla^{s} u+\frac{\lambda}{2\mu}\left({\partial_{t}\sigma}+ u\cdot\nabla\sigma\right)\nonumber\end{align} (3.3) \begin{align} -\frac{\lambda}{2\mu}\left(\sigma\cdot\nabla u+\left(\nabla u\right)^{T}\cdot\sigma\right) &={0}\qquad\textrm{in}\:\varOmega,\:t\in\left]0,T\right[,\end{align} (3.3) \begin{align} u &={0}\qquad\textrm{on}\:\partial\varOmega,\:t\in\left]0,T\right[,\end{align} (3.4) \begin{align} u\mid_{t=0}\ = u_{0}\qquad\!\!\textrm{in}\:\varOmega, \end{align} (3.5) \begin{align} \sigma\mid_{t=0}\ =\sigma_{0}\qquad\!\!\textrm{in}\:\varOmega. \end{align} (3.6) The unknowns of the problem are the velocity field $$u\left ({x},t\right )$$, the pressure $$p\left ({x},t\right )$$ and the viscoleastic or elastic part $$\sigma \left ({x},t\right )$$ of the extra stress tensor. The physical parameters are the dynamic viscosity $$\mu$$, the density of the fluid $$\rho$$, a real parameter $$\beta \in \left [0,1\right ]$$ to define the amount of viscous or solvent viscosity $$\left (\mu _{s}=\beta \mu \right )$$ and elastic or polymeric viscosity $$\left (\mu _{p}=(1-\beta )\mu \right )$$ in the fluid, and the relaxation time $$\lambda$$ that represents the elasticity of the fluid. Finally, $$f\in (H^{-1}(\varOmega ))^{d}$$ is the external volume force applied to the fluid confined in $$\varOmega$$. For viscoelastic fluids, the problem is incomplete without the definition of a constitutive equation for the elastic stresses $$\left (\sigma \right )$$. A large variety of approaches exist to define it (see Bird et al., 1987a,b for a complete description). In this work, we use the classical Oldroyd-B constitutive model (3.3) for this purpose. The conservation laws (3.1)–(3.2) together with the Oldroyd-B constitutive equation (3.3) are a mixed parabolic-hyperbolic problem that needs to be complemented with initial (3.5)–(3.6) and boundary (3.4) conditions to close the problem. For simplicity in the exposition, we will consider the simplest boundary condition for the velocity and no boundary conditions for the stress field. With respect to the initial conditions, $$u_{0}$$ and $$\sigma _{0}$$ are functions defined on the whole domain $$\varOmega$$. For a complete description of the mathematical structure of the problem, we refer to Fernández-Cara et al. (2002) and Renardy (1989). 3.2 Variational form To write the weak form of problem (3.1)–(3.3) we need to introduce some functional spaces. Let $${\mathscr{V}}=({H_{0}^{1}}(\varOmega ))^{d}$$, $${\varUpsilon }:=\left \{ \tau \mid \tau \in (L^{2}\left (\varOmega \right ))_{\mathrm{sym}}^{d\times d},~{w}\cdot \nabla \tau \in (L^{2}(\varOmega ))^{d\times d}_{\mathrm{sym}} ~\forall \,{w}\in{\mathscr{V}}\right \}$$ (subscript sym standing for symmetric tensors) and $$Q=L^{2}(\varOmega )/\mathbb{R}$$, the spaces of the velocity, the elastic stresses and the pressure, respectively. If we denote $$U:=\left (u,p,\sigma \right )$$, $${\mathscr{X}}:={\mathscr{V}}\times Q\times{\varUpsilon }$$, the weak form of the problem consists in finding $$U:\left ]0,T\right [\rightarrow{\mathscr{X}}$$ such that the initial conditions are satisfied and \begin{align} \rho\left({\partial_{t} u}, v\right)+\rho\langle u\cdot\nabla u, v\rangle+\left(2\beta\mu\nabla^{s} u+\sigma,\nabla^{s} v\right)-\left(p,\nabla\cdot v\right) =\left\langle f, v\right\rangle\!, \end{align} (3.7) \begin{align}\qquad\qquad\qquad\qquad\qquad \left(\nabla\cdot u,q\right) =0, \end{align} (3.8) \begin{align} \left(\frac{1}{2\mu}\sigma-(1-\beta)\nabla^{s} u,\tau\right)+\frac{\lambda}{2\mu}\left({\partial_{t}\sigma}+ u\cdot\nabla\sigma,\tau\right)-\frac{\lambda}{2\mu}\left(g\left(u,\sigma\right),\tau\right) =0, \end{align} (3.9) for all $${V}=\left (v,q,\tau \right )\in{\mathscr{X}}$$. The last term of the constitutive equation represents the traction or rotational term, defined as $$g\left (u,\sigma \right ):=\sigma \cdot \nabla u+\left (\nabla u\right )^{T}\cdot \sigma$$. In a compact form, problem (3.7)–(3.9) can be written as: \begin{align} \left(\mathscr{D}_{t}\left(U\right),{V}\right)+B\left(u,\sigma;U,{V}\right)=\left\langle f, v\right\rangle\!, \end{align} (3.10) for all $${V}\in{\mathscr{X}}$$, where \begin{align} B\left(\hat{u},\hat{\sigma};U,{V}\right)= &\, \rho\langle\hat{u}\cdot\nabla u, v\rangle+2\beta\mu\left(\nabla^{s} u,\nabla^{s} v\right)+\left(\sigma,\nabla^{s} v\right)-\left(p,\nabla\cdot v\right)+\left(\nabla\cdot u,q\right)\\ & +\frac{1}{2\mu}\left(\sigma,\tau\right)-(1-\beta)\left(\nabla^{s} u,\tau\right)+\frac{\lambda}{2\mu}\left(\hat{u}\cdot\nabla\sigma,\tau\right)-\frac{\lambda}{2\mu}\left(g\left(\hat{u},\hat{\sigma}\right),\tau\right)\!,\nonumber \end{align} (3.11) and \begin{align} \left(\mathscr{D}_{t}\left(U\right),{V}\right):=\rho\left({\partial_{t} u}, v\right)+\frac{\lambda}{2\mu}\left({\partial_{t}\sigma},\tau\right)\!. \end{align} (3.12) 3.3 Galerkin finite element discretization From $$\mathscr{T}_{h}$$ we may construct conforming finite element spaces for the velocity, the pressure and the elastic stress, $${\mathscr{V}}_{h}\subset{\mathscr{V}}$$, $$Q_{h}\subset Q$$ and $${\varUpsilon }_{h}\subset{\varUpsilon }$$, respectively. Denoting $${\mathscr{X}}_{h}={\mathscr{V}}_{h}\times Q_{h}\times{\varUpsilon }_{h}$$, the Galerkin finite element approximation of problem (3.10) consists in finding $$U_{h}:\left ]0,T\right [\rightarrow{\mathscr{X}}_{h}$$ such that \begin{align} \left(\mathscr{D}_{t}\left(U_{h}\right),{V}_{h}\right)+B\left(u_{h},\sigma_{h};U_{h},{V}_{h}\right)=\left\langle \,f, v_{h}\right\rangle\!, \end{align} (3.13) for all $${V}_{h}\in{\mathscr{X}}_{h}$$, and satisfying the appropriate initial conditions. Until now, we have posed no restrictions on the choice of the finite element spaces. However, let us analyse the numerical stability of problem (3.13). If we take $${V}_{h}=U_{h1}=\left ((1-\beta ) u_{h},(1-\beta )p_{h},\sigma _{h}\right )$$, it is found that \begin{align} B\left({u}_{h},{\sigma}_{h};U_{h},U_{h1}\right) =2(1-\beta)\beta\mu\left\Vert \nabla^{s} u_{h}\right\Vert{}^{2}+\frac{1}{2\mu}\left\Vert \sigma_{h}\right\Vert{}^{2} -\frac{\lambda}{2\mu}\left(g(u_{h},\sigma_{h}),\sigma_{h}\right)\!. \end{align} (3.14) It is seen from (3.14) that B is not coercive in $${\mathscr{X}}_{h}$$, and we can ensure only control on $$\left \Vert \sigma _{h}\right \Vert$$ for all $$\beta$$ assuming $$\lambda \nabla u_{h}$$ to be small enough. To ensure the control of $$p_{h}$$ and $$\nabla ^{s} u_{h}$$, one has then to choose finite element spaces satisfying \begin{align} \underset{q_{h}\in{{Q}}_{h}}{\textrm{inf}}\,\underset{v_{h}\in{\mathscr{V}}_{h}}{\textrm{sup}}\frac{\left(q_{h},\nabla\cdot v_{h}\right)}{\left\Vert v_{h}\right\Vert{}_{{\mathscr{V}}_{h}}\left\Vert q_{h}\right\Vert{}_{{{Q}}_{h}}} \geqslant C_{1}, \end{align} (3.15) \begin{align} \underset{v_{h}\in{\mathscr{V}}_{h}}{\textrm{inf}}\,\underset{\tau_{h}\in\mathscr{{\varUpsilon}}_{h}}{\textrm{sup}}\frac{\left(\tau_{h},\nabla^{s} v_{h}\right)}{\left\Vert \tau_{h}\right\Vert{}_{\mathscr{{\varUpsilon}}_{h}}\left\Vert v_{h}\right\Vert{}_{{\mathscr{V}}_{h}}} \geqslant C_{2}, \end{align} (3.16) where $$C_{1}$$ and $$C_{2}$$ are positive constants. It is therefore required that the finite element spaces satisfy (3.15) and (3.16), which is a stringent requirement inherited from the mixed form of the Navier–Stokes problem (Castillo & Codina, 2014a). The two compatibility conditions of the viscoelastic flow problem do not allow us the use of an arbitrary interpolation for the different fields because the scheme may become unstable. The implementation of inf–sup stable elements is a possible solution for this problem; however, from the numerical point of view, the spaces that fulfil these conditions are limited and complex, particularly when the problem needs to be solved in three dimensions. The other possibility is to use a stabilized formulation that permits the use of any interpolation for the variables, which is the approach studied in this work. Note that the constitutive equation is of convective nature, and therefore, some kind of stabilization technique has to be used even if inf–sup stable elements are used, and likewise for the momentum equation. 4. Stabilized finite element method In general, a stabilized formulation consists of replacing B in (3.10) by another multilinear form $$B_{\textrm{stab}}$$, possibly mesh dependent, designed to enhance stability without upsetting accuracy. The formulation we analyse, proposed in the study by Castillo & Codina (2014b), is described next. For the sake of conciseness, let us consider equal order continuous interpolation for all variables. The case of discontinuous pressures and stresses can be treated using the technique employed in the study by Castillo & Codina (2017b). Thus, let us consider $${\mathscr{X}}_{h}={\mathscr{V}}_{h}\times Q_{h}\times{\varUpsilon }_{h}$$, where $${\mathscr{V}}_{h}=\left [V_{h}\cap{H_{0}^{1}}(\varOmega )\right ]^{d}$$, $$Q_{h}=[V_{h}\cap \mathscr{C}^{0}(\varOmega )]/\mathbb{R}$$ and $${\varUpsilon }_{h}=\left [V_{h}\cap \mathscr{C}^{0}(\varOmega )\right ]_{\textrm{sym}}^{d\times d}$$. The method consists in replacing (3.13) by the following problem: find $$U_{h} :\ ]0, T[\ \rightarrow{\mathscr{X}}_{h}$$ such that \begin{align} \left(\mathscr{D}_{t}\left(U_{h}\right),{V}_{h}\right)+B_{\textrm{stab}}\left(u_{h},\sigma_{h};U_{h},{V}_{h}\right)=\left\langle \,f, v\right\rangle\!, \end{align} (4.1) for all $${V}_{h}\in{\mathscr{X}}_{h}$$, where \begin{align} B_{\textrm{stab}}\left(\hat{u}_{h},\hat{\sigma}_{h};U_{h},{V}_{h}\right)=B\left(\hat{u}_{h},\hat{\sigma}_{h};U_{h},{V}_{h}\right)+B^{*}\left(\hat{u}_{h},\hat{\sigma}_{h};U_{h},{V}_{h}\right)\!, \end{align} (4.2) and $$B^{*}$$represents the additional stabilization terms added to the Galerkin formulation. Using the same notation as in the study by Castillo & Codina (2014b), we can define $$B^{*}$$ as $$B^{*}\left(\hat{u}_{h},\hat{\sigma}_{h};U_{h},{V}_{h}\right)= S_{1}^{\perp}\left(\hat{u}_{h};U_{h},{V}_{h}\right)+S_{2}^{\perp}\left(U_{h},{V}_{h}\right)+S_{3}^{\perp}\left(\hat{u}_{h},\hat{\sigma}_{h};U_{h},{V}_{h}\right)\!,$$ where \begin{align} S_{1}^{\perp}\left(\hat{u}_{h};U_{h},{V}_{h}\right)= &\, \alpha_{u}(\hat{u}_{h})\left(P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla u_{h}\right],P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla v_{h}\right]\right)+\alpha_{u}(\hat{u}_{h})\left(P_{u}^{\perp}\left[\nabla p_{h}\right],P_{u}^{\perp}\left[\nabla q_{h}\right]\right)\nonumber\\ & +\alpha_{u}(\hat{u}_{h})\left(P_{u}^{\perp}\left[\nabla\cdot\sigma_{h}\right],P_{u}^{\perp}\left[(1-\beta)\nabla\cdot\tau_{h}\right]\right)\!, \end{align} (4.3) \begin{align} S_{2}^{\perp}\left(U_{h},{V}_{h}\right)= \alpha_{p}(\hat{u}_{h})\left(P_{p}^{\perp}\left[\nabla\cdot u_{h}\right],P_{p}^{\perp}\left[\nabla\cdot v_{h}\right]\right)\!, \end{align} (4.4) \begin{align} S_{3}^{\perp}\left(\hat{u}_{h},\hat{\sigma}_{h};U_{h},{V}_{h}\right)=\, & \alpha_{\sigma}(\hat{u}_{h})\left(P_{\sigma}^{\perp}\left[-(1-\beta)\nabla^{s} u_{h}+\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla\sigma_{h}-g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right)\right]\right.\!,\nonumber\\ & \left.P_{\sigma}^{\perp}\left[-\nabla^{s} v_{h}+\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla\tau_{h}+g^{*}\left(\hat{u}_{h},\tau_{h}\right)\right)\right]\right)\!, \end{align} (4.5) and $$g^{*}\left (\hat{u}_{h},\tau _{h}\right )$$ represents the adjoint operator of $$g\left (\hat{u}_{h},\tau _{h}\right )$$, defined as $$g^{*}\left (\hat{u}_{h},\tau _{h}\right ):=\tau _{h}\cdot \left (\nabla \hat{u}_{h}\right )^{T}+ \ \nabla \hat{u}_{h}\cdot \tau _{h}$$. Here $$P_{u}^{\perp }=I-P_{u}$$, where $$P_{u}:L^{2}(\varOmega )\rightarrow \tilde{{\mathscr{V}}}_{h}$$ is the $$L^{2}(\varOmega )$$-projection onto $$\tilde{{\mathscr{V}}}_{h}$$, the velocity space without boundary conditions, and $$P_{p}^{\perp }$$ and $$P_{\sigma }^{\perp }$$ are defined in an analogous way. We will also need $$P_{u,0}:=L^{2}(\varOmega )\rightarrow{{\mathscr{V}}}_{h}$$, that is to say, the $$L^{2}(\varOmega )$$-projection onto the velocity space with boundary conditions. Note that all the stabilization terms in (4.3)–(4.5) are multiplied by $$\alpha _{i}$$, $$i=u,p,\sigma$$. These terms are the components of the stabilization parameter matrix, that can be defined as $${\alpha}=\operatorname{diag}\left(\alpha_{u}{I}_{d},\alpha_{p},\alpha_{\sigma}{I}_{d\times d}\right)\!,$$ where \begin{align} \alpha_{u}(\hat{u}_{h})= \left[c_{1}\frac{\mu}{h^{2}}+c_{2}\frac{\rho\left\Vert \hat{u}_{h}\right\Vert{}_{\infty}}{h}\right]^{-1}, \end{align} (4.6) \begin{align} \alpha_{p}(\hat{u}_{h})= \frac{h^{2}}{c_{1}\alpha_{u}(\hat{u}_{h})}, \end{align} (4.7) \begin{align} \alpha_{\sigma}(\hat{u}_{h})= \left[c_{3}\frac{1}{2\mu}+\frac{\lambda}{2\mu}\left(c_{4}\frac{\left\Vert \hat{u}_{h}\right\Vert{}_{\infty}}{h}+2c_{5}\left\Vert \nabla\hat{u}_{h}\right\Vert{}_{\infty}\right)\right]^{-1}, \end{align} (4.8) where $$c_{i}$$, i = 1, 2, 3, 4, 5, are algorithmic constants and $${I}_{d}$$ and $${I}_{d\times d}$$ are the second and fourth order identity tensors, respectively. A general approach to design the terms of the stabilization parameter matrix was proposed in the study by Codina (2009) for the three-field Stokes problem. In this work, it is shown that the parameters can be uniquely determined by dimensionality, assuming that this matrix is diagonal. In general, to get optimal control the stabilization parameters need to be evaluated elementwise (or even pointwise). The constant expression adopted in this work allows us to simplify the analysis. The stabilizing mechanisms introduced by the terms $$S_{1}^{\perp }$$, $$S_{2}^{\perp }$$ and $$S_{3}^{\perp }$$ are the following. The first component of $$S_{1}^{\perp }$$ gives control on the convective term, the second component gives control on the pressure gradient and the third term gives control on the divergence of the viscoelastic stress. The term $$S_{2}^{\perp }$$ is not a must, but in some cases it improves stability of the problem. Finally, the term $$S_{3}^{\perp }$$ adds stability in the constitutive equation. Note that some of the components of this last term are the convective–convective term of the viscoelastic stress tensor and an equivalent EVSS-structure component, amongst other cross local inner product terms (see Castillo & Codina, 2014b for more details of this spatial stabilized formulation). The addition of these three terms permits the approximation of convection dominated problems both in velocity and in stress and the implementation of equal order interpolation for all the unknowns. The orthogonal projections introduce consistency errors, but of optimal order, a key point in the design of accurate nonresidual-based methods. For stationary problems, the resulting formulation turns out to have optimal order of convergence, as checked numerically in the study by Castillo & Codina (2014b) for linear and quadratic elements. The term-by-term form of $$S_{1}^{\perp }$$ was proposed instead of a residual-based one because the former shows a better numerical behaviour in problems where high gradients in pressure and stress are present (see Castillo & Codina, 2017a for more details about this fact). We will need a condition on the interpolating spaces that holds in the case of equal order interpolations (see Codina & Blasco, 2000) and that can be written as follows: given $${a}_{h}, v_{h}\in{{\mathscr{V}}}_{h}$$, $$q_{h}\in{{{Q}}}_{h}$$, $$\tau _{h}\in{\varUpsilon }_{h}$$ and $${z}_{h}:=\rho{a}_{h}\cdot \nabla v_{h}+\nabla q_{h}-\nabla \cdot \tau _{h}$$, there holds \begin{align} \left\Vert{z}_{h}\right\Vert \leqslant c_{\mathrm{m}}\left(\left\Vert P_{u,0}\left({z}_{h}\right)\right\Vert +\left\Vert P_{u}^{\perp}\left({z}_{h}\right)\right\Vert \right), \end{align} (4.9) for a constant $$c_{\mathrm{m}}>0$$. According to this condition, the component of $$P_{u}\left ({z}_{h}\right )$$ that corresponds to the boundary of $$\varOmega$$ can be bounded in terms of the right-hand side (RHS) of equation (4.9). To prove this, one can use the macroelement technique employed in the study by Codina & Blasco (2000). 5. Linearized time-dependent case The numerical analysis of the stabilized formulation presented in this work is divided in two steps. In this section we present the stability analysis of the linearized case. The second part (Section 6), is devoted to the nonlinear analysis. 5.1 Linearized stabilized semidiscrete problem The equations for incompressible viscoelastic flows have several nonlinear terms, both in the momentum and in the constitutive equation. In the former we have the convective term and in the latter we have the term corresponding to the convection of stresses and the rotational term arising from the objective derivative of stresses. On top of that, the stabilization terms depend also on the velocity, introducing therefore additional nonlinearities. As it is usual for incompressible flow problems, for the convective term of the momentum equation we can use a Picard scheme as linearization technique, taking the advection velocity as the velocity of the previous iteration. This leads only to first order convergence, but it is a robust option. The constitutive equation is rather more complex, and sometimes the implementation of combined algorithms is recommended. See for example the work of Castillo & Codina (2014b), where a Newton scheme was combined with a continuation method, or the work of Howell (2009), where different types of continuation methods were proposed. For simplicity in the numerical analysis, we will use only the fixed point scheme for all the nonlinear terms, including the terms associated to the matrix of stabilization parameters $${\alpha }$$. Therefore, we analyse the following semidiscrete linearized problem: given $$\hat{u}_{h}:\ ]0,T[\ \rightarrow{\mathscr{V}}_{h}$$ and $$\hat{\sigma }_{h}: \ ]0,T[\ \rightarrow{\varUpsilon }_{h}$$, find $$U_{h}:\left ]0,T\right [\rightarrow{\mathscr{X}}_{h}$$ such that \begin{align} \rho\left({\partial_{t} u_{h}},v_{h}\right)&+\rho\tilde{c}\left(\hat{u}_{h}, u_{h}, v_{h}\right)+\left(2\beta\mu\nabla^{s}u_{h}+\sigma_{h},\nabla^{s} v_{h}\right)-\left(p_{h},\nabla\cdot v_{h}\right)\nonumber\\ &+\alpha_{u}\left(P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla u_{h}\right],P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla v_{h}\right]\right)+\alpha_{p}\left(P_{p}^{\perp}\left[\nabla\cdot u_{h}\right],P_{p}^{\perp}\left[\nabla\cdot v_{h}\right]\right)\nonumber \\ &+\alpha_{\sigma}\left(P_{\sigma}^{\perp}\left[(1-\beta)\nabla^{s}u_{h}-\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla\sigma_{h}-g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right)\right],P_{\sigma}^{\perp}\left[\nabla^{s} v_{h}\right]\right) =\left\langle\, f, v_{h}\right\rangle\!, \end{align} (5.1) \begin{align} \qquad\qquad\qquad\qquad\qquad\qquad\;\left(\nabla\cdot u_{h},q_{h}\right)+\alpha_{u}\left(P_{u}^{\perp}\left[\nabla p_{h}\right],P_{u}^{\perp}\left[\nabla q_{h}\right]\right) =0, \end{align} (5.2) \begin{align} \left(\frac{1}{2\mu}\sigma_{h}-(1-\beta)\nabla^{s} u_{h},\tau_{h}\right)&+\frac{\lambda}{2\mu}\left({\partial_{t}\sigma_{h}},\tau_{h}\right)+\frac{\lambda}{2\mu}\tilde{c}\left(\hat{u}_{h},\sigma_{h},\tau_{h}\right)\nonumber\\ &-\frac{\lambda}{2\mu}\left(g\left(\hat{u}_{h},\hat{\sigma}_{h}\right),\tau_{h}\right)+(1-\beta)\alpha_{u}\left(P_{u}^{\perp}\left[\nabla\cdot\sigma_{h}\right],P_{u}^{\perp}\left[\nabla\cdot\tau_{h}\right]\right)\nonumber \\ &+\alpha_{\sigma}\left(P_{\sigma}^{\perp}\left[-(1-\beta)\nabla^{s}u_{h}+\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla\sigma_{h}-g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right)\right],\right.\nonumber \\ &\quad\left.P_{\sigma}^{\perp}\left[\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla\tau_{h}+g^{*}\left(\hat{u}_{h},\tau_{h}\right)\right)\right]\right) =0, \end{align} (5.3) for all $${V}_{h}=\left (v_{h},q_{h},\tau _{h}\right )\in{\mathscr{X}}_{h}$$. The initial conditions are set as appropriate projections of $$u_{0}$$ and $$\sigma _{0}$$. The stabilization parameters are computed using $$\hat{u}_{h}$$ (see (4.6)–(4.8)). 5.2 Existence and uniqueness of the semidiscrete solution The following existence and uniqueness analysis of the discrete solution was motivated by the procedure followed in the study by Burman & Fernández (2007) for the two-field Navier–Stokes problem. To prove existence and uniqueness of the discrete linearized problem (5.1)–(5.3), we shall make use of the following pressure and velocity subspaces: \begin{align*} Q_{h}^{\star}= & \left\{ q_{h}\in Q_{h}\vert\,\left(P_{u,0}^{\perp}\left[\nabla q_{h}\right],P_{u,0}^{\perp}\left[\nabla q_{h}\right]\right)=0\right\}\!,\\{\mathscr{V}}_{h}^{\operatorname{div}}= & \left\{ v_{h}\in{\mathscr{V}}_{h}\vert\,\left(q_{h},\nabla\cdot v_{h}\right)=0,\:\forall\, q_{h}\in Q_{h}^{\star}\right\}\!. \end{align*} In addition, $$Q_{h}\backslash Q_{h}^{\star }$$ will stand for the supplementary of $$Q_{h}^{\star }$$ in $$Q_{h}$$, i.e., $$Q_{h}=\left (Q_{h}\backslash Q_{h}^{\star }\right )\oplus Q_{h}^{\star }$$. To ensure that $${\mathscr{V}}_{h}^{\operatorname{div}}$$ is not trivial, we use the following lemma: Lemma 5.1 There exists a constant $$\gamma>0$$, independent of h, such that $$\underset{q_{h}\in{{Q}}_{h}^{\star}}{\textrm{inf}}\,\underset{v_{h}\in{\mathscr{V}}_{h}}{\textrm{sup}}\frac{\left(q_{h},\nabla\cdot v_{h}\right)}{\left\Vert v_{h}\right\Vert{}_{1}\left\Vert q_{h}\right\Vert}\geqslant\gamma\!.$$ Proof. Let $$q_{h}\in Q_{h}^{\star }$$. From inf–sup theory (see for example Girault & Raviart, 1986), there exists $$v_{q}\in ({H_{0}^{1}}(\varOmega ))^{d}$$ such that $$\nabla\cdot v_{q}=q_{h},\;\left\Vert v_{q}\right\Vert{}_{1}\lesssim\left\Vert q_{h}\right\Vert\!.$$ Integrating by parts, we have \begin{align*} \left\Vert q_{h}\right\Vert{}^{2}= & \left(q_{h},\nabla\cdot v_{q}\right)\\ = & -\left(\nabla q_{h}, v_{q}-P_{u,0}\left[v_{q}\right]\right)+\left(q_{h},\nabla\cdot P_{u,0}\left[v_{q}\right]\right)\\ = & \left(q_{h},\nabla\cdot P_{u,0}\left[v_{q}\right]\right)\!, \end{align*} since, as $$q_{h}\in{{Q}}_{h}^{\star }$$, then $$\nabla q_{h}\in{\mathscr{V}}_{h}$$, and $$v_{q}-P_{u,0}\left [v_{q}\right ]$$ belongs to the orthogonal to this space. In addition, using the quasi-uniformity of the mesh and the Poincaré–Friedrics inequality, we have $$\left\Vert P_{u,0}\left[v_{q}\right]\right\Vert{}_{1}^{2}\lesssim\left\Vert \nabla v_{q}\right\Vert{}^{2}\lesssim\left\Vert q_{h}\right\Vert{}^{2},$$ which completes the proof. Theorem 5.2 (Existence-uniqueness) The semidiscrete problem (5.1)–(5.3) has a unique solution. Proof. Problem (5.1)–(5.3), satisfying the initial conditions, can be written in operator form as \begin{align}{3} M_{u}\partial_{t} u_{h}+K_{u}\left(\hat{u}_{h}\right) u_{h}+Gp_{h}-D_{\sigma}\sigma_{h} =M_{u} f\; \quad \textrm{in}\:{\mathscr{V}}_{h}, \end{align} (5.4) \begin{align} D u_{h} =S_{1,p}^{\perp}p_{h}\;\; \textrm{in}\:Q_{h}, \end{align} (5.5) \begin{align} M_{\sigma}\partial_{t}\sigma_{h}+K_{\sigma}\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\sigma_{h}-S u_{h} =0\qquad\;\; \textrm{in}\:{\varUpsilon}_{h}. \end{align} (5.6) These are equations posed in terms of linear forms defined on the spaces specified in each equation. The components of these forms when applied to $$v_{h}\in{\mathscr{V}}_{h}$$, $${q}_{h}\in{{{Q}}}_{h}$$ and $$\tau _{h}\in{{\varUpsilon }}_{h}$$ are given by \begin{align*} M_{u} f (v_{h})= &\, \langle f, v_{h}\rangle,\\ (K_{u}(\hat{u}_{h}) u_{h}) (v_{h}) = &\, \rho\tilde{c}\left(\hat{u}_{h}, u_{h}, v_{h}\right)+2\beta\mu\left(\nabla^{s} u_{h}\nabla^{s} v_{h}\right)+\alpha_{u}\left(P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla u_{h}\right],P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla v_{h}\right]\right)\\ & +\alpha_{p}\left(P_{p}^{\perp}\left[\nabla\cdot u_{h}\right],P_{p}^{\perp}\left[\nabla\cdot v_{h}\right]\right)\\ & +(1-\beta)\alpha_{\sigma}\left(P_{\sigma}^{\perp}\left[(1-\beta)\nabla^{s} u_{h}-\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla\sigma_{h}-g\left(\hat{u}_{h},\sigma_{h}\right)\right)\right],P_{\sigma}^{\perp}\left[\nabla^{s} v_{h}\right]\right),\\ Gp_{h}(v_{h})= & -\left(p_{h},\nabla\cdot v_{h}\right)\!,\\ D_{\sigma}\sigma_{h}(v_{h})= & \left(\sigma_{h},\nabla^{s}v_{h}\right)\!,\\ D u_{h} (q_{h}) = & \left(\nabla\cdot u_{h},q_{h}\right)\!,\\ S_{1,p}^{\perp}p_{h} (q_{h}) = & -\alpha_{u}\left(P_{u}^{\perp}\left[\nabla q_{h}\right],P_{u}^{\perp}\left[\nabla q_{h}\right]\right)\!,\\ M_{\sigma}\partial_{t}\sigma_{h} (\tau_{h})= &\, \frac{\lambda}{2\mu}\left({\partial_{t}\sigma_{h}},\tau_{h}\right)\!,\\ (K_{\sigma}(\hat{u}_{h},\hat{\sigma}_{h})\sigma_{h})(\tau_{h})= &\, \frac{1}{2\mu}\left(\sigma_{h},\tau_{h}\right)+\frac{\lambda}{2\mu}\tilde{c}\left(\hat{u}_{h},\sigma_{h},\tau_{h}\right)+(1-\beta)\alpha_{u}\left(P_{u}^{\perp}\left[\nabla\cdot\sigma_{h}\right],P_{u}^{\perp}\left[\nabla\cdot\tau_{h}\right]\right)\\ & +\alpha_{\sigma}\left(P_{\sigma}^{\perp}\left[-(1-\beta)\nabla^{s}u_{h}+\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla\sigma_{h}-g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right)\right],\right.\\ &\quad \left.P_{\sigma}^{\perp}\left[\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla\tau_{h}+g^{*}\left(\hat{u}_{h},\tau_{h}\right)\right)\right]\right)\!,\\ S u_{h}(\tau_{h})= & -(1-\beta)\left(\nabla^{s} u_{h},\tau_{h}\right)\!. \end{align*} We also introduce the operator $$D^{1}:{\mathscr{V}}_{h} \longrightarrow ({Q}_{h}^{\star })^{\prime }$$, defined by $$D^{1}v_{h}= \left(Dv_{h}\right)_{Q_{h}^{\star}},\:\:\:\forall\, v_{h}\in{\mathscr{V}}_{h},$$ where the subscript stands for the restriction to $${Q_{h}^{\star }}$$. From Lemma 5.1, it follows that $$D^{1}$$ is surjective and $$\left (D^{1}\right )^{T}$$ is injective, and therefore, $${\mathscr{V}}_{h}^{\operatorname{div}}:=\ker \left (D^{1}\right )\neq \left \{\textbf{0}\right \}$$. Let us consider the following reduced formulation, derived from (5.1)–(5.3) with $${V}_{h}\in{\mathscr{X}}_{h}^{\star }={\mathscr{V}}_{h}^{\operatorname{div}}\times \left (Q_{h}\backslash Q_{h}^{\star }\right )\times{\varUpsilon }_{h}$$: find $$U_{h} = (u_{h},\tilde{p}_{h},\sigma _{h}) : \left ]0,T\right [\rightarrow{\mathscr{X}}_{h}^{\star }$$ such that \begin{align} M_{u}\partial_{t} u_{h}+K_{u}\left(\hat{u}_{h}\right) u_{h}+G\tilde{p}_{h}-D_{\sigma}\sigma_{h} =M_{u} f\qquad \textrm{in}\:{\mathscr{V}}_{h}^{\operatorname{div}}, \end{align} (5.7) \begin{align} D u_{h} =S_{1,p}^{\perp}\tilde{p}_{h}\quad\; \textrm{in}\:\left(Q_{h}\backslash Q_{h}^{\star}\right)\!, \end{align} (5.8) \begin{align} M_{\sigma}\partial_{t}\sigma_{h}+K_{\sigma}\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\sigma_{h}-S u_{h} =0\qquad\;\;\;\; \textrm{in}\:{\varUpsilon}_{h}. \end{align} (5.9) Since, by construction, $$Q_{h}^{\star }=\operatorname{Ker}\left (S_{1,p}^{\perp }\right )$$, we conclude that $$S_{1,p}^{\perp }$$ is invertible in $$Q_{h}\backslash Q_{h}^{\star }$$, and therefore, from (5.8) we obtain $$\tilde{p}_{h}=\left(S_{1,p}^{\perp}\right)_{Q_{h}\backslash Q_{h}^{\star}}^{-1}D u_{h}.$$ By plugging this expression into (5.7), we obtain the following equivalent problem: \begin{align} M_{u}\partial_{t} u_{h}+K_{u}\left(\hat{u}_{h}\right) u_{h}+G\left(S_{1,p}^{\perp}\right)_{Q_{h}\backslash Q_{h}^{\star}}^{-1}D u_{h}-D_{\sigma}\sigma_{h} =M_{u} f\:\:\: \textrm{in}\:{\mathscr{V}}_{h}^{\operatorname{div}}, \end{align} (5.10) \begin{align} M_{\sigma}\partial_{t}\sigma_{h}+K_{\sigma}\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\sigma_{h}-S u_{h} =0\quad\;\;\; \textrm{in}\:{\varUpsilon}_{h}, \end{align} (5.11) which is a standard Cauchy problem for $$u_{h}$$ and $$\sigma _{h}$$. Existence and uniqueness for $$u_{h}$$ and $$\sigma _{h}$$ follows by the Lipschitz continuity of $$K_{u}\left (\hat{u}_{h}\right )$$ and $$K_{\sigma }\left (\hat{u}_{h},\hat{\sigma }_{h}\right )$$. We may then recover $$\tilde{p}_{h}$$ uniquely from (5.8). Therefore, the reduced problem (5.7)–(5.9) has a unique solution. On the other hand, from (5.7) it follows that $$M_{u}\partial_{t} u_{h}+K_{u}\left(\hat{u}_{h}\right)u_{h}+G\tilde{p}_{h}-D_{\sigma}\sigma_{h}-M_{u} f \in\left(\operatorname{Ker}\left(D^{1}\right)\right)^{0},$$ with $$\left (\operatorname{Ker}\left (D^{1}\right )\right )^{0}$$ standing for the polar set of $$\operatorname{Ker}\left (D^{1}\right )$$. From Lemma 5.1, it follows that $$D^{1}$$ is an isomorphism from $$Q_{h}^{\star }$$ onto $$\left (\operatorname{Ker}\left (D^{1}\right )\right )^{0}$$. Thus, there exists a unique $$p^{1}\in Q_{h}^{\star }$$ such that \begin{align} M_{u}\partial_{t}u_{h}+K_{u}\left(\hat{u}_{h}\right) u_{h}+G\tilde{p}_{h}-D_{\sigma}\sigma_{h}-M_{u} f= Gp^{1}\:\:\:\textrm{in}\,{\mathscr{V}}_{h}. \end{align} (5.12) Therefore, from (5.12) and the reduced problem (5.7)–(5.9), it follows that the discrete problem (5.1)–(5.3) has a unique solution, given by $$\left (u_{h},\tilde{p}_{h}-p^{1},\sigma _{h}\right )$$. 5.3 Stability of the linearized semidiscrete problem We will use the following working norm in the subsequent analysis: \begin{align} \left\Vert{V}_{h}\right\Vert{}_{W}^{2}= \, & 2\beta\mu(1-\beta)\left\Vert \nabla^{s}v_{h}\right\Vert{}^{2}+\frac{1}{2\mu}\left\Vert \tau_{h}\right\Vert{}^{2}+(1-\beta)\alpha_{u}\left\Vert P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla v_{h}\right]\right\Vert^{2}\nonumber \\ & +(1-\beta)\alpha_{u}\left\Vert \rho{\partial_{t} v_{h}}+\rho\hat{u}_{h}\cdot\nabla v_{h}+\nabla q_{h}-\nabla\cdot\tau_{h}\right\Vert^{2}\nonumber\\ & +(1-\beta)\alpha_{p}\left\Vert \nabla\cdot v_{h}\right\Vert{}^{2}+(1-\beta)\alpha_{u}\left\Vert P_{u}^{\perp}\left[\nabla q_{h}\right]\right\Vert^{2}+(1-\beta)\alpha_{u}\left\Vert P_{u}^{\perp}\left[\nabla\cdot\tau_{h}\right]\right\Vert^{2}\nonumber \\ & +\alpha_{\sigma}\left\Vert \frac{\lambda}{2\mu}{\partial_{t}\tau_{h}}-(1-\beta)\nabla^{s}v_{h}+\frac{\lambda}{2\mu}\hat{u}_{h}\cdot\nabla\tau_{h}\right\Vert^{2}. \end{align} (5.13) The term multiplied by $$\alpha _{p}$$ is not strictly necessary, since it is already contained in $$\left \Vert \nabla ^{s} v_{h}\right \Vert$$, but sometimes could reinforce stability. We will keep it for generality, to see the effect of the stabilizing term associated to the pressure. From the definition of $$\Vert \cdot \Vert _{W}$$ it is observed that we have some control on the convective terms of the momentum equation and of the constitutive equation. In view of the expression of the stabilization parameters, this control remains meaningful in the convection dominated limit for the momentum equation. This property will no longer be valid in the nonlinear case, since, as it is standard in the analysis of nonlinear problems involving the Navier–Stokes equation, the results will be proved in the diffusion dominated regime. The next result states the stability of the proposed semidiscrete formulation, defined in (5.1)–(5.3). Theorem 5.3 (Stability). Let $$\left (u_{h},p_{h},\sigma _{h}\right )$$ be the solution of (5.1)–(5.3). Then, the following stability holds, for almost all $$t\in \left [0,T\right ]$$: \begin{align} & (1-\beta)\frac{\rho}{2}\left\Vert u_{h}\left(t\right)\right\Vert{}^{2}+\frac{\lambda}{4\mu}\left\Vert \sigma_{h}\left(t\right)\right\Vert{}^{2} +{\int_{0}^{t}}\left\Vert (u_{h},p_{h},\sigma_{h}) \right\Vert{}_{W}^{2}\,\mathrm{d}t \nonumber \\ & \qquad \lesssim \frac{c_{\mathrm{K}}^{2}}{2\mu}{\int_{0}^{t}}\left\Vert f\right\Vert{}_{H^{-1}}^{2}\,\mathrm{d}t+\frac{\lambda^{2}}{2\mu}{\int_{0}^{t}}\left\Vert g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right\Vert{}^{2}\,\mathrm{d}t +(1-\beta)\frac{\rho}{2}\left\Vert u_{0}\right\Vert{}^{2}+\frac{\lambda}{4\mu}\left\Vert \sigma_{0}\right\Vert{}^{2}. \end{align} (5.14) Proof. In this proof $$\varepsilon _{1},\varepsilon _{2},...$$ are positive constants used in the application of Young’s inequalities. The values will be chosen at the end of the proof. Taking $$\left (v_{h},q_{h},\tau _{h}\right )=U_{h1}=\left ((1-\beta )u_{h},(1-\beta )p_{h},\sigma _{h}\right )$$ in (5.1)–(5.3), adding up the resulting equations and using Cauchy–Schwarz’s and Young’s inequalities, we arrive at \begin{align} & (1-\beta)\frac{\rho}{2}\frac{\mathrm{d}}{\mathrm{d}t}\left\Vert u_{h}\right\Vert{}^{2}+2\mu\beta(1-\beta)\left\Vert \nabla^{s} u_{h}\right\Vert{}^{2}+\frac{\lambda}{4\mu}\frac{\mathrm{d}}{\mathrm{d}t}\left\Vert \sigma_{h}\right\Vert{}^{2}\nonumber\\ &\quad +(1-\beta)\alpha_{p}\left\Vert P_{p}^{\perp}\left[\nabla\cdot u_{h}\right]\right\Vert^{2}+(1-\beta)\alpha_{u}\left\Vert P_{u}^{\perp}\left[\nabla\cdot\sigma_{h}\right]\right\Vert^{2}\nonumber \\ & \quad+\frac{1}{2\mu}\left(1-\frac{\varepsilon_{1}}{2}-4\frac{\alpha_{\sigma}}{2\mu}\lambda^{2}\left\Vert \nabla\hat{u}_{h}\right\Vert{}_{\infty}^{2}\left(\frac{1}{2\varepsilon_{2}}+\frac{1}{2\varepsilon_{4}}\right)\right)\left\Vert \sigma_{h}\right\Vert{}^{2}\nonumber \\ & \quad+(1-\beta)\alpha_{u}\left\Vert P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla u_{h}\right]\right\Vert^{2}+(1-\beta)\alpha_{u}\left\Vert P_{u}^{\perp}\left[\nabla p_{h}\right]\right\Vert^{2}\nonumber \\ & \quad+\alpha_{\sigma}\left(1-\frac{\varepsilon_{2}}{2}-\frac{\varepsilon_{3}}{2}\right)\left\Vert P_{\sigma}^{\perp}\left[-(1-\beta)\nabla^{s} u_{h}+\frac{\lambda}{2\mu}\left({\partial_{t}\sigma_{h}}+\hat{u}_{h}\cdot\nabla\sigma_{h}\right)\right]\right\Vert^{2}\nonumber \\ &\qquad \leqslant (1-\beta)\left\langle f, u_{h}\right\rangle +\frac{\lambda^{2}}{4\mu}\left[\frac{1}{\varepsilon_{1}}\left\Vert g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right\Vert{}^{2}+\frac{\alpha_{\sigma}}{\mu}\left(\frac{1}{2\varepsilon_{3}}+\frac{\varepsilon_{4}}{2}\right)\left\Vert P_{\sigma}^{\perp}\left[g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right]\right\Vert^{2}\right]. \end{align} (5.15) If we compare the bounded terms of expression (5.15) with the terms of the working norm (5.13), we can see that the missing terms are all of them associated to the finite element space. The key point is that this missing control comes from the Galerkin part of the multilinear form $$B_{\operatorname{stab}}$$ in equation (4.2). Let us take $${V}_{h1}=\left (v_{h},q_{h},\tau _{h}\right )=\alpha _{u}\left ((1-\beta )v_{1},0,{0}\right )$$ with $$v_{1}\equiv P_{u,0}\left[\rho{\partial_{t} u_{h}}+\rho\hat{u}_{h}\cdot\nabla u_{h}+\nabla p_{h}-\nabla\cdot\sigma_{h}\right]\!.$$ Then, we can proceed in a similar way as in the $$U_{h1}$$ case. Taking $$v_{1}$$ as test function, integrating by parts the terms arising from the divergence of the total stress and the gradient of the pressure, using the fact that $$\hat{u}_{h}={0}$$ on $$\partial \varOmega$$, and applying Cauchy–Schwarz’s and Young’s inequalities and inverse inequalities, we obtain the following result \begin{align} &(1-\beta)\alpha_{u}\left(1 -\alpha_{u}c_{\operatorname{inv}}^{2}\left[\frac{\alpha_{p}}{h^{2}}\frac{\varepsilon_{6}}{2}+\frac{\rho\left\Vert \hat{u}_{h}\right\Vert{}_{\infty}}{h}\frac{\varepsilon_{7}}{2}+\frac{\alpha_{\sigma}}{h^{2}}\left(\frac{\varepsilon_{8}}{2}+\frac{\varepsilon_{9}}{2}\right)\right]\right.\nonumber\\ &\quad \left. -\beta\alpha_{u}\frac{2\mu}{h^{2}}c_{\operatorname{inv}}^{2}\frac{\varepsilon_{5}}{2}-\alpha_{u}\frac{d}{4\varepsilon_{10}}\frac{c_{\mathrm{K}}^{2}}{2\mu}\rho^{2}\left\Vert \nabla\hat{u}_{h}\right\Vert{}_{\infty}^{2}\right)\left\Vert v_{1}\right\Vert{}^{2}\nonumber \\ &\quad -2\mu(1-\beta)\left(\beta\frac{1}{2\varepsilon_{5}}+\frac{\varepsilon_{10}}{2}\right)\left\Vert \nabla^{s} u_{h}\right\Vert{}^{2}-(1-\beta)\alpha_{p}\frac{1}{2\varepsilon_{6}}\left\Vert P_{p}^{\perp}\left[\nabla\cdot u_{h}\right]\right\Vert^{2}\nonumber \\ &\quad -(1-\beta){\alpha_{u}^{2}}\frac{\rho\left\Vert \hat{u}_{h}\right\Vert{}_{\infty}}{h}\frac{1}{2\varepsilon_{7}}\left\Vert P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla u_{h}\right]\right\Vert^{2}\nonumber \\ &\quad -(1-\beta)\alpha_{\sigma}\frac{1}{2\varepsilon_{8}}\left\Vert P_{\sigma}^{\perp}\left[\frac{\lambda}{2\mu}{\partial_{t}\sigma_{h}}-(1-\beta)\nabla^{s} u_{h}+\frac{\lambda}{2\mu}\hat{u}_{h}\cdot\nabla\sigma_{h}\right]\right\Vert^{2}\nonumber \\ &\quad\quad\leqslant (1-\beta)\alpha_{u}\left\langle f, v_{1}\right\rangle +\frac{(1-\beta)}{2\varepsilon_{9}}\frac{\lambda^{2}}{2\mu}\frac{\alpha_{\sigma}}{2\mu}\left\Vert P_{\sigma}^{\perp}\left[g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right]\right\Vert^{2}. \end{align} (5.16) To obtain the control of $$P_{p}\left [\nabla \cdot u_{h}\right ]$$, we proceed taking $${V}_{h2}=(1-\beta )\alpha _{p}\left ({0},q_{2},{0}\right )$$ with $$q_{2}\equiv P_{p}\left [\nabla \cdot u_{h}\right ]$$: \begin{align} (1-\beta)\alpha_{p}\left(\left(\nabla\cdot u_{h},q_{2}\right)+\alpha_{u}\left(P_{u}^{\perp}\left[\nabla p_{h}\right],\nabla q_{2}\right)\right)\geqslant &\, (1-\beta)\alpha_{p}\left(1-c_{\operatorname{inv}}^{2}\alpha_{u}\frac{\alpha_{p}}{h^{2}}\frac{\varepsilon_{11}}{2}\right)\left\Vert P_{p}\left[\nabla\cdot u_{h}\right]\right\Vert{}^{2}\nonumber \\ & -(1-\beta)\alpha_{u}\frac{1}{2\varepsilon_{11}}\left\Vert P_{u}^{\perp}\left[\nabla p_{h}\right]\right\Vert^{2}. \end{align} (5.17) Finally, taking: $${V}_{h3}=\alpha _{\sigma }\left ({0},0,\sigma _{3}\right )$$ with $$\sigma _{3}\equiv P_{\sigma }\left [\frac{\lambda }{2\mu }{\partial _{t}\sigma _{h}}-(1-\beta )\nabla ^{s} u_{h}+\frac{\lambda }{2\mu }\hat{u}_{h}\cdot \nabla \sigma _{h}\right ]$$, and proceeding as before, we obtain \begin{align} & -\frac{1}{2\mu}\left(\frac{1}{2\varepsilon_{12}}+\frac{1}{2}\frac{1}{2\varepsilon_{17}}\right)\left\Vert \sigma_{h}\right\Vert{}^{2}-\frac{(1-\beta)}{2\varepsilon_{14}}\alpha_{u}\left\Vert P_{u}^{\perp}\left[\nabla\cdot\sigma_{h}\right]\right\Vert^{2}\nonumber\\ & +\alpha_{\sigma}\left(1-\frac{\alpha_{\sigma}}{2\mu}\frac{\varepsilon_{12}}{2}-\frac{\varepsilon_{13}}{2}-\frac{\varepsilon_{17}}{4}\frac{\alpha_{\sigma}}{2\mu}d\left(\lambda\left\Vert \nabla\hat{u}_{h}\right\Vert{}_{\infty}\right)^{2}-c_{\operatorname{inv}}^{2}(1-\beta)\alpha_{u}\frac{\alpha_{\sigma}}{h^{2}}\frac{\varepsilon_{14}}{2}\vphantom{\left[\left(\frac{\lambda\left\Vert \hat{u}_{h}\right\Vert{}_{\infty}}{2\mu h}\right)^{2}+4\left(\frac{\lambda\left\Vert \nabla\hat{u}_{h}\right\Vert{}_{\infty}}{2\mu}\right)^{2}\right]}\right.\nonumber \\ & \left.-\,c_{\operatorname{inv}}^{2}2\alpha_{\sigma}^{2}\left[\left(\frac{\lambda\left\Vert \hat{u}_{h}\right\Vert{}_{\infty}}{2\mu h}\right)^{2}+4\left(\frac{\lambda\left\Vert \nabla\hat{u}_{h}\right\Vert{}_{\infty}}{2\mu}\right)^{2}\right]\left(\frac{1}{2\varepsilon_{15}}+\frac{1}{2\varepsilon_{16}}\right)\right)\left\Vert \sigma_{3}\right\Vert{}^{2}\nonumber \\ & -\alpha_{\sigma}\frac{\varepsilon_{15}}{2}\left\Vert P_{\sigma}^{\perp}\left[\frac{\lambda}{2\mu}{\partial_{t}\sigma_{h}}-(1-\beta)\nabla^{s} u_{h}+\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla\sigma_{h}\right)\right]\right\Vert^{2}\nonumber \\ & \quad \leqslant \frac{\lambda^{2}}{4\mu}\frac{\alpha_{\sigma}}{\mu}\left(\frac{1}{2\varepsilon_{13}}\left\Vert g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right\Vert{}^{2}+\frac{\varepsilon_{16}}{2}\left\Vert P_{\sigma}^{\perp}\left[g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right]\right\Vert^{2}\right). \end{align} (5.18) Let $${V}_{h}^{*}=U_{h1}+\theta _{1}{V}_{h1}+\theta _{2}{V}_{h2}+\theta _{3}{V}_{h3}$$, with $${V}_{hi}$$, i = 1, …, 3 introduced above. Adding (5.16)–(5.18) multiplied by $$\theta _{i}$$, i = 1, …, 3, and adding also (5.15), we arrive at an expresion of the form $$LHS\left({V}_{h}^{*}\right)\leqslant RHS\left({V}_{h}^{*}\right)\!.$$ For the RHS, applying Young’s inequalities and the inverse inequality, we obtain \begin{align} RHS\left({V}_{h}^{*}\right) &\leqslant (1-\beta)\left(\frac{1}{2\mu}\left(\frac{1}{2\varepsilon_{17}}+\theta_{1}\frac{1}{2\varepsilon_{18}}\right)c_{\varOmega}^{2}\left\Vert f\right\Vert{}_{H^{-1}}^{2}\right.\nonumber\\ & \quad\left.+\ 2\mu\frac{\varepsilon_{17}}{2}\left\Vert \nabla^{s} u_{h}\right\Vert{}^{2}+\theta_{1}\alpha_{u}\left(\alpha_{u}\frac{2\mu}{h^{2}}\right)c_{\operatorname{inv}}^{2}\frac{\varepsilon_{18}}{2}\left\Vert v_{1}\right\Vert{}^{2}\right)\nonumber \\ & \quad+\frac{\lambda^{2}}{2\mu}\left(\frac{1}{2\varepsilon_{1}}+\theta_{3}\frac{\alpha_{\sigma}}{2\mu}+\frac{\alpha_{\sigma}}{2\mu}\left(\frac{1}{2\varepsilon_{3}}+\frac{\varepsilon_{4}}{2}+\theta_{1}\frac{(1-\beta)}{2\varepsilon_{9}}+\theta_{3}\frac{\varepsilon_{16}}{2}\right)\right)\left\Vert g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right\Vert{}^{2}. \end{align} (5.19) For the left-hand side (LHS), integrating inequalities (5.15), (5.16), (5.17), (5.18) and (5.19) from 0 to t, we obtain \begin{align} & (1-\beta)\frac{\rho}{2}\left\Vert u_{h}\right\Vert{}^{2}\left(t\right)+\frac{\lambda}{4\mu}\left\Vert \sigma_{h}\right\Vert{}^{2}\left(t\right)+2\mu(1-\beta)\beta C_{1}{\int_{0}^{t}}\left\Vert \nabla^{s}u_{h}\right\Vert{}^{2}\,\mathrm{d}t\nonumber\\ &\quad +\frac{1}{2\mu}C_{2}{\int_{0}^{t}}\left\Vert \sigma_{h}\right\Vert{}^{2}\,\mathrm{d}t+(1-\beta)C_{3}\alpha_{u}{\int_{0}^{t}}\left\Vert P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla u_{h}\right]\right\Vert^{2}\,\mathrm{d}t\nonumber \\ &\quad +(1-\beta)C_{4}\alpha_{p}{\int_{0}^{t}}\left\Vert P_{p}^{\perp}\left[\nabla\cdot u_{h}\right]\right\Vert^{2}\,\mathrm{d}t+(1-\beta)C_{5}\alpha_{u}{\int_{0}^{t}}\left\Vert P_{u}^{\perp}\left[\nabla\cdot\sigma_{h}\right]\right\Vert^{2}\,\mathrm{d}t\nonumber \\ &\quad +(1-\beta)C_{6}\alpha_{u}{\int_{0}^{t}}\left\Vert P_{u}^{\perp}\left[\nabla p_{h}\right]\right\Vert^{2}\,\mathrm{d}t+C_{7}\alpha_{\sigma}{\int_{0}^{t}}\left\Vert P_{\sigma}^{\perp}\left[-(1-\beta)\nabla^{s} u_{h}+\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla\sigma_{h}\right)\right]\right\Vert^{2}\,\mathrm{d}t\nonumber \\ &\quad +(1-\beta)C_{8}\alpha_{u}{\int_{0}^{t}}\left\Vert P_{u,0}\left[\rho{\partial_{t} u_{h}}+\rho\hat{u}_{h}\cdot\nabla u_{h}+\nabla p_{h}-\nabla\cdot\sigma_{h}\right]\right\Vert{}^{2}\,\mathrm{d}t\nonumber \\ &\quad +(1-\beta)C_{9}\alpha_{p}{\int_{0}^{t}}\left\Vert P_{p}\left[\nabla\cdot u_{h}\right]\right\Vert{}^{2}\,\mathrm{d}t\nonumber \\ &\quad +C_{10}\alpha_{\sigma}{\int_{0}^{t}}\left\Vert P_{\sigma}\left[\frac{\lambda}{2\mu}{\partial_{t}\sigma_{h}}-(1-\beta)\nabla^{s}u_{h}+\frac{\lambda}{2\mu}\hat{u}_{h}\cdot\nabla\sigma_{h}\right]\right\Vert^{2}\,\mathrm{d}t\nonumber \\ & \qquad \leqslant (1\!-\!\beta)\frac{c_{\varOmega}^{2}}{2\mu}C_{11}{\int_{0}^{t}}\!\left\Vert f\right\Vert{}_{H^{-1}}^{2}\,\mathrm{d}t\!+\!\frac{\lambda^{2}}{2\mu}C_{12}{\int_{0}^{t}}\!\left\Vert g\!\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\!\right\Vert{}^{2}\,\mathrm{d}t+(1-\beta)\frac{\rho}{2}\left\Vert u_{h}\left(0\right)\right\Vert{}^{2}\!+\!\frac{\lambda}{4\mu}\!\left\Vert \sigma_{h}\left(0\right)\right\Vert{}^{2}\!, \end{align} (5.20) where $$C_{i}$$, i = 1, …, 12, can be easily identified. The result then follows by choosing $$\varepsilon _{1},\ldots ,\varepsilon _{18}$$ and $$\theta _{1}, \theta _{2} ,\theta _{3}$$, in such a way that $$C_{i}>0$$ for all i, which is possible by choosing some of the constants and the algorithmic stabilization parameters sufficiently small. The missing control on $$\rho \partial _{t}u_{h}+\rho \hat{u}_{h}\cdot \nabla u_{h}+\nabla p_{h}-\nabla \cdot \sigma _{h}$$ follows from (4.9). 6. Analysis of the nonlinear problem In this section we show that under suitable conditions, a solution to the discretized system exists. Using the same procedure proposed by Ervin & Miles (2003), the proof can be subdivided in the following four steps: Define an iterative map in such a way that a fixed point of the map is a solution to the original problem. Show that the map is well defined and bounded on bounded sets. Show that there exists an invariant ball of the map. Apply Schauder’s fixed point theorem to conclude there exists a discrete solution in this ball. The results to be obtained in this section yield existence of a semidiscrete solution, as well as stability and convergence. Due to the hypotheses needed in the proof of the results (most of which are already present in the analysis of the Navier–Stokes equation), the norm in which the results are shown is weaker than the norm used in the linearized problem. In essence, we will have $$L^{\infty }(0,T)$$-control for both the $$L^{2}(\varOmega )$$-norm of velocity and stresses and $$L^{2}(0,T)$$-control for the $$H^{1}(\varOmega )$$-norm of velocity. The pressure, on the other hand, is controlled only in a norm involving the stabilization term, and not the natural $$L^{2}(\varOmega )$$-norm. The precise assumptions required on the continuous solution are: Assumption 6.1 System (3.10) has a solution $$\left (u,p,\sigma \right )$$ continuous in time and satisfying \begin{alignat}{3} & \underset{0\leqslant t\leqslant T}{\sup}\left\Vert u\right\Vert{}_{\infty}\leqslant D_{1}, \quad & \underset{0\leqslant t\leqslant T}{\sup}\left\Vert \nabla u\right\Vert{}_{\infty}\leqslant D_{2}, \quad & \underset{0\leqslant t\leqslant T}{\sup}\left\Vert u\right\Vert{}_{k+1}\leqslant D_{3},\nonumber\\ & \underset{0\leqslant t\leqslant T}{\sup}\left\Vert \sigma\right\Vert{}_{\infty}\leqslant D_{4}, \quad & \underset{0\leqslant t\leqslant T}{\sup}\left\Vert \nabla\sigma\right\Vert{}_{\infty}\leqslant D_{5}, \quad & \underset{0\leqslant t\leqslant T}{\sup}\left\Vert \sigma\right\Vert{}_{k+1}\leqslant D_{6}, \\ & \underset{0\leqslant t\leqslant T}{\sup}\left\Vert p\right\Vert{}_{k}\leqslant D_{7}, \quad & \underset{0\leqslant t\leqslant T}{\sup}\left\Vert \partial_{t} u\right\Vert{}_{k}\leqslant D_{8}, \quad & \underset{0\leqslant t\leqslant T}{\sup}\left\Vert \partial_{t}\sigma\right\Vert{}_{k}\leqslant D_{9},\nonumber \end{alignat} (6.1) for certain positive constants $$D_{i}$$i = 1, …, 9, which are supposed to be small enough. For the time-discrete problem, if $$\delta t$$ is the time step size, one usually needs a condition of the form $$\delta t \geqslant C \alpha _{u}$$ for a positive constant C, which is encountered in most stabilized finite element methods; see Codina et al. (2007) and Badia & Codina (2009) and references therein for a description of the problem and a way to avoid this restriction, which we shall not consider in this work. In the time continuous case, the boundedness in time of $$\Vert p\Vert _{k}$$ and the assumption that T is large enough allows us to prove convergence, as we show next. We do not pretend however to consider the long-term behaviour of the solution, which would require the modification of the stabilized formulation and the analysis presented in the study by Badia et al. (2010). In order to write all estimates in dimensionless form, let $$L_{d}$$ be a characteristic length of the problem and $$T_{d}$$ a characteristic time scale. These parameters may explode with the viscosity, and therefore the estimates are not valid for high Reynolds numbers. The main result on existence and convergence reads as follows: Theorem 6.2 (Convergence) Let k be the interpolation order, assumed to be the same for all variables, with $$k\geqslant d/2$$. Suppose also that Assumption 6.1 holds, that T is sufficiently large and that the $$L^{2}(0,T; (H^{-1}(\varOmega ))^{d})$$-norm of f is bounded. Then, if the viscosity is sufficiently large, there exists a solution to (5.1)–(5.3) satisfying \begin{alignat}{1} \underset{0\leqslant t\leqslant T}{\sup}\left\Vert u-u_{h}\right\Vert{}^{2}+\frac{1}{T_{d}}{\int_{0}^{T}}\left(L_{d}\left\Vert \nabla\left(u-u_{h}\right)\right\Vert \right)^{2}\,\mathrm{d}t\leqslant \, &\, u_{\star}^{2}h^{2k},\nonumber\\ \underset{0\leqslant t\leqslant T}{\sup}\left\Vert \sigma-\sigma_{h}\right\Vert{}^{2}\leqslant \, &\, \sigma_{\star}^{2}h^{2k}, \\{\int_{0}^{T}}\alpha_{u} \left\Vert \nabla\left(p-p_{h}\right)\right\Vert{}^{2}\,\mathrm{d} t\leqslant \, &\, p_{\star}^{2}h^{2k}, \nonumber \end{alignat} (6.2) where $$u_{\ast }$$, $$p_{\ast }$$ and $$\sigma _{\ast }$$ are appropriate dimensional factors that render the estimates dimensionally consistent. Proof. Step 1. The iterative map. A mapping $$\delta :L^{2}\left (0,T;{\mathscr{V}}_{h}\right )\times L^{2}\left (0,T;Q_{h}\right )\times L^{\infty }\left (0,T;{\varUpsilon }_{h}\right )\rightarrow L^{2}\left (0,T;{\mathscr{V}}_{h}\right )\times L^{2}\left (0,T;Q_{h}\right )\times L^{\infty }\left (0,T;{\varUpsilon }_{h}\right )$$ is defined via $$\left (u_{h},p_{h},\sigma _{h}\right )=\delta \left (\hat{u}_{h},\hat{p}_{h},\hat{\sigma }_{h}\right )$$, where $$U_{h} = \left (u_{h},p_{h},\sigma _{h}\right )$$ satisfies (5.1)–(5.3), for all $${V}_{h}=\left (v_{h},q_{h},\tau _{h}\right )\in{\mathscr{X}}_{h}$$. Thus, given an initial guess for the three unknowns $$\hat{U}_{h} := \left (\hat{u}_{h},\hat{p}_{h},\hat{\sigma }_{h}\right )$$, solving the above system for $$\left (u_{h},p_{h},\sigma _{h}\right )$$ gives a new approximation to the solution. Also, it is clear that the fixed point is a solution to the approximating system (5.1)–(5.3), i.e., $$\delta \left ({u}_{h},{p}_{h},{\sigma }_{h}\right )=\left ({u}_{h},{p}_{h},{\sigma }_{h}\right )$$ implies that $$\left ({u}_{h},{p}_{h},{\sigma }_{h}\right )$$ is a solution to (4.1). Step 2. The mapping $$\delta$$ is well defined and bounded on bounded sets. The existence and uniqueness of the discrete solution was proved in subsection 5.2 for a linearized problem, that can be associated to the solution of the fixed-point problem. The stability result proved in subsection 5.3 ensures that the linearized problem, that can be viewed as a fixed-point iteration of the nonlinear case, is stable and bounded under suitable regularity assumptions. Note that this step is crucial in the definition of the fixed point mapping and can be used to establish that the mapping $$\delta$$ is bounded on bounded sets. Step 3. Existence of an invariant ball. We begin defining an invariant ball. Let $$R=h^{k}$$, and for $$V = (v,q,\tau )$$ let us define the norm \begin{align} \Vert V {\Vert_{B}^{2}} := \rho \sup_{0\leqslant t\leqslant T}\Vert v\Vert^{2} + \frac{\rho}{T_{d}}{\int_{0}^{T}} (L_{d}\Vert \nabla v \Vert)^{2}\,\mathrm{d} t + \frac{\lambda}{\mu} \sup_{0\leqslant t\leqslant T} \Vert \tau\Vert^{2} + \alpha_{u} {\int_{0}^{T}} \Vert \nabla q \Vert^{2}\,\mathrm{d} t. \end{align} (6.3) As we shall see, our result could be proved in a finer norm, but the additional terms that we could include in $$\Vert \cdot \Vert _{B}$$ (see below) do not provide significant additional control for large viscosities. Let us define now the ball $$\mathscr{B}_{h}$$ as \begin{align} \mathscr{B}_{h}= \Bigl\{ V_{h} = \left(v_{h},q_{h},\tau_{h}\right) : \ ]0,T[\ \rightarrow{\mathscr{X}}_{h}~~\textrm{such that}~~ \Vert V_{h} -U\Vert_{B} \leqslant U_{\ast} R \Bigr\}, \end{align} (6.4) where $$U^{2}_{\ast } = \rho u^{2}_{\ast } + \frac{\lambda }{\mu }\sigma ^{2}_{\ast } + p^{2}_{\ast }$$ and $$u_{\ast }$$, $$p_{\ast }$$ and $$\sigma _{\ast }$$ are constructed along the proof. In the definition (6.4), $$U = \left (u,p,\sigma \right )$$ is the solution of (3.10). Let us pick now $$\hat{U}_{h} := \left (\hat{u}_{h},\hat{p}_{h},\hat{\sigma }_{h}\right )\in \mathscr{B}_{h}$$, arbitrary, and let $$U_{h} = \left (u_{h},p_{h},\sigma _{h}\right )=\delta \left (\hat{u}_{h},\hat{p}_{h},\hat{\sigma }_{h}\right )$$; it satisfies (5.1)–(5.3). Then, it is readily checked that \begin{align} \rho\left(\partial_{t}\left(u-u_{h}\right), v_{h}\right)&+\rho\left(\tilde{c}\left(u,u, v_{h}\right)-\tilde{c}\left(\hat{u}_{h},u_{h}, v_{h}\right)\right)\nonumber\\ &+2\beta\mu\left(\nabla^{s}\left(u-u_{h}\right),\nabla^{s} v_{h}\right)+\left(\sigma-\sigma_{h},\nabla^{s} v_{h}\right)-\left(p-p_{h},\nabla\cdot v_{h}\right)\nonumber \\ &-\alpha_{u}\left(P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla u_{h}\right],P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla v_{h}\right]\right)-\alpha_{p}\left(P_{p}^{\perp}\left[\nabla\cdot u_{h}\right],P_{p}^{\perp}\left[\nabla\cdot v_{h}\right]\right)\nonumber \\ &+\left(\nabla\cdot\left(u-u_{h}\right),q_{h}\right)-\alpha_{u}\left(P_{u}^{\perp}\left[\nabla p_{h}\right],P_{u}^{\perp}\left[\nabla q_{h}\right]\right)\nonumber \\ &+\frac{1}{2\mu}\left(\sigma-\sigma_{h},\tau_{h}\right)-(1-\beta)\left(\nabla^{s}\left(u-u_{h}\right),\tau_{h}\right)+\frac{\lambda}{2\mu}\left(\partial_{t}\left(\sigma-\sigma_{h}\right),\tau_{h}\right)\nonumber \\ &+\frac{\lambda}{2\mu}\left(\tilde{c}\left(u,\sigma,\tau_{h}\right)-\tilde{c}\left(\hat{u}_{h},\sigma_{h},\tau_{h}\right)\right)-\frac{\lambda}{2\mu}\left(g\left(u,\sigma\right),\tau_{h}\right)+\frac{\lambda}{2\mu}\left(g\left(\hat{u}_{h},\hat{\sigma}_{h}\right),\tau_{h}\right)\nonumber \\ &-(1-\beta)\alpha_{u}\left(P_{u}^{\perp}\left[\nabla\cdot\sigma_{h}\right],P_{u}^{\perp}\left[\nabla\cdot\tau_{h}\right]\right)\nonumber \\ &-\alpha_{\sigma}\left(P_{\sigma}^{\perp}\left[-(1-\beta)\nabla^{s}u_{h}+\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla\sigma_{h}-g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right)\right],\right.\nonumber \\ &\quad\left.P_{\sigma}^{\perp}\left[-\nabla^{s} v_{h} + \frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla\tau_{h}+g^{*}\left(\hat{u}_{h},\tau_{h}\right)\right)\right]\right) =0. \end{align} (6.5) Let us define the following approximation and interpolation errors: \begin{align*} \varLambda=u-{\mathscr{U}} & \:\:\textrm{and}\:\: E={\mathscr{U}}-u_{h},\\ \varGamma=\sigma-{\varSigma} & \:\:\textrm{and}\:\: F={\varSigma}-\sigma_{h},\\ \varPi=p-\mathscr{P} & \:\:\textrm{and}\:\:G=\mathscr{P}-p_{h}, \end{align*} where $${\mathscr{U}}=P_{u,0}\left [u\right ]$$, $${\varSigma }=P_{\sigma }\left [\sigma \right ]$$ and $$\mathscr{P}=P_{p}\left [p\right ]$$. We obviously have that $${e}_{u} = u-u_{h}=\varLambda +E$$, $${e}_{\sigma } = \sigma -\sigma _{h}=\varGamma + F$$, $$e_{p} = p-p_{h}=\varPi +G$$. Let us describe the strategy to prove that U belongs to the ball $$\mathscr{B}_{h}$$. Let $$P[U] = ({\mathscr{U}},{\varSigma }, \mathscr{P})$$. From (6.5) we will show that \begin{align} \Vert P[U] - U_{h}{\Vert_{B}^{2}} \leqslant \varphi_{1}(D) \Vert U -\hat{U}_{h}{\Vert^{2}_{B}} + \Vert U - P[U]{\Vert^{2}_{I}}, \end{align} (6.6) where $$\varphi _{1}(D)$$ is a certain function polynomial in terms of the components of the array of constants $$D = (D_{1},\dots ,D_{9})$$ introduced in Assumption 6.1 and $$\Vert \cdot \Vert _{I}$$ is a norm of what we may consider the interpolation error U − P[U]; this norm will appear along the proof. Thanks to the properties of the interpolation, we will check that \begin{align} \Vert U - P[U]{\Vert_{I}^{2}} \leqslant \varphi_{2}(D) h^{2k}, \end{align} (6.7) and we can already verify that \begin{align} \Vert U - P[U]{\Vert_{B}^{2}} \leqslant \varphi_{3}(D) h^{2k}. \nonumber \end{align} Again, $$\varphi _{2}(D)$$ and $$\varphi _{3}(D)$$ are functions polynomial in terms of the components of D. Using the triangle inequality, (6.6)–(6.7) and the fact that $$\hat{U}_{h} \in \mathscr{B}_{h}$$, we will have that \begin{align} \Vert U - U_{h}{\Vert_{B}^{2}} & \leqslant \Vert P[U] - U_{h}{\Vert_{B}^{2}} + \Vert U - P[U]{\Vert_{B}^{2}} \nonumber\\ & \leqslant \varphi_{1}(D) U^{2}_{\ast} h^{2k} + \varphi_{2}(D) h^{2k} + \varphi_{3}(D) h^{2k}.\nonumber \end{align} Taking the components of array D sufficiently small, we will be able to guarantee that \begin{align} \Vert U - U_{h}{\Vert_{B}^{2}} \leqslant U^{2}_{\ast} h^{2k},\nonumber \end{align} that is to say, $${U}_{h} \in \mathscr{B}_{h}$$. Therefore, the goal is to prove (6.6) and check that (6.7) holds. The constant $$U_{\ast }$$, from which we can obtain $$u_{\ast }$$, $$p_{\ast }$$ and $$\sigma _{\ast }$$, can be taken as $$U^{2}_{\ast } = c( \varphi _{2}(D) + \varphi _{3}(D))$$, with c > 1, provided D is such that $$\varphi _{1}(D) + c^{-1} \leqslant 1$$. Taking $$v_{h}=(1-\beta )\left (E+\gamma \alpha _{u}{P_{u,0}\left [\nabla G\right ]}\right )$$, with $$\gamma> 0$$, $$\tau _{h}= F$$ and $$q_{h}=(1-\beta )G$$ in (6.5), we obtain \begin{align} \rho(1-\beta) \frac{\mathrm{d}}{\mathrm{d}t}\left\Vert E\right\Vert{}^{2}+\frac{\lambda}{\mu} \frac{\mathrm{d}}{\mathrm{d}t}\left\Vert F\right\Vert{}^{2}+ Q(t) \lesssim N(t) + S(t) +{\overset{23}{\underset{j=1}{\sum}}R_{j}}(t), \end{align} (6.8) where \begin{align*} Q(t) & := \beta(1-\beta)\mu\left\Vert E\right\Vert{}_{1}^{2} +\frac{1}{\mu}\left\Vert F\right\Vert{}^{2}+(1-\beta)\alpha_{u}\left\Vert P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla E\right]\right\Vert^{2}+(1-\beta)\alpha_{p}\left\Vert P_{p}^{\perp}\left[\nabla\cdot E\right]\right\Vert^{2}\nonumber \\ &\qquad +(1-\beta)\alpha_{u}\left(\left\Vert P_{u}^{\perp}\left[\nabla G\right]\right\Vert{}^{2}+{\gamma\left\Vert{P_{u,0}\left[\nabla G\right]}\right\Vert{}^{2}}\right)+(1-\beta)\alpha_{u}\left\Vert P_{u}^{\perp}\left[\nabla\cdot F\right]\right\Vert^{2}\nonumber \\ & \qquad + \alpha_{\sigma}\Bigl\Vert P_{\sigma}^{\perp} \Bigl[ -(1-\beta)\nabla^{s} E+\frac{\lambda}{2\mu}\hat{u}_{h}\cdot\nabla F\Bigr] \Bigr\Vert^{2},\nonumber\\ N(t) & := -\rho(1-\beta)\left(\tilde{c}\left(u,u, E\right)-\tilde{c}\left(\hat{u}_{h},u_{h}, E\right)\right)-\frac{\lambda}{2\mu}\left(\tilde{c}\left(u,\sigma,F\right)-\tilde{c}\left(\hat{u}_{h},\sigma_{h},F\right)\right)\nonumber \\ &\qquad {-\,(1-\beta)\gamma\alpha_{u}\rho\left(\tilde{c}\left(u,u,{P_{u,0}\left[\nabla G\right]}\right)+\tilde{c}\left(\hat{u}_{h},u_{h},{P_{u,0}\left[\nabla G\right]}\right)\right)},\nonumber \\ S(t) & := -\,\alpha_{\sigma}\Bigl(P_{\sigma}^{\perp}\Bigl[-(1-\beta)\nabla^{s} E+\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla F\right)\Bigr],g^{*}(\hat{u}_{h},F)\Bigr), \end{align*} and the terms $$R_{j}(t)$$ are defined below. From (6.8), and using (4.9) for the pressure term, we have that \begin{align} \Vert P[U] - U_{h}{\Vert^{2}_{B}} & \lesssim \rho(1-\beta) \sup_{0\leqslant t\leqslant T} \Vert E\Vert^{2} + \frac{\lambda}{\mu} \sup_{0\leqslant t\leqslant T} \Vert F\Vert^{2} + {\int_{0}^{T}} Q(t) \,\mathrm{d}t \nonumber\\ &\lesssim{\int_{0}^{T}} \Biggl(N(t) + S(t) +{\overset{23}{\underset{j=1}{\sum}}R_{j}}(t) \Biggr)\,\mathrm{d}t. \end{align} (6.9) The objective is to see that all these terms in the RHS can be bounded as indicated in (6.6) or absorbed by the LHS. We will not detail the positive dimensionless constants that appear in this process, in which we will make frequent use of Young’s inequality; the parameter that appears when using it will be generically denoted by $$\varepsilon$$ or $$\xi$$, understanding that it is small enough. We could track the different appearances of these parameters as we have done in the proof of Theorem 5.3, and choose them at the end; however, since we are not interested in the values of the constants, we will proceed in a more conceptual way. We will start with the terms $$R_{j}(t)$$ of (6.8). Using the $$L^{2}$$-orthogonality between $$\varLambda$$ and E and between $$\varGamma$$ and F and Young’s inequality, we have that \begin{align*} R_{1}(t) & =-(1-\beta)\rho\left(\partial_{t}\varLambda, E\right) = 0,\\ R{}_{2}(t) & =-\frac{\lambda}{\mu}\left(\partial_{t}\varGamma,F\right) = 0,\\ R_{3}(t) & =(1-\beta)\left(\varPi,\nabla\cdot E\right)\lesssim \mu(1-\beta){\varepsilon}\left\Vert E\right\Vert{}_{1}^{2}+\frac{1}{\mu}\frac{(1-\beta)}{\varepsilon}\left\Vert \varPi\right\Vert{}^{2}. \end{align*} The term that involves the discrete error of the pressure needs a special treatment. Using (4.9) we have that \begin{align*} R_{4}(t) & =-(1-\beta)\left(\nabla\cdot\varLambda,G\right)= (1-\beta)\left(\varLambda,\nabla G\right)\\ & \lesssim \alpha_{u}^{-1}\frac{(1-\beta)}{\varepsilon}\left\Vert \varLambda\right\Vert{}^{2}+(1-\beta){\varepsilon}\alpha_{u}\left( \left\Vert{P_{u}^{\perp}\left[\nabla G\right]}\right\Vert{}^{2} + \left\Vert{{P_{u,0}\left[\nabla G\right]}}\right\Vert{}^{2} \right). \end{align*} The following terms only need Cauchy–Schwarz’s and Young’s inequalities: \begin{align*} R_{5}(t) & =-2\beta(1-\beta)\mu\left(\nabla^{s}\varLambda,\nabla^{s} E\right)\lesssim \beta(1-\beta)\mu{\varepsilon}\left\Vert E\right\Vert{}_{1}^{2}+ \beta(1-\beta)\mu\frac{1}{\varepsilon}\left\Vert \varLambda\right\Vert{}_{1}^{2},\\ R_{6}(t) & =-(1-\beta)\left(\varGamma,\nabla^{s} E\right)\lesssim \mu(1-\beta){\varepsilon}\left\Vert E\right\Vert{}_{1}^{2}+\frac{1}{\mu}(1-\beta)\frac{1}{\varepsilon}\left\Vert \varGamma\right\Vert{}^{2},\\ R_{7}(t) & =-\frac{1}{2\mu}\left(\varGamma,F\right) = 0,\\ R_{8}(t) & =(1-\beta)\left(\nabla^{s}\varLambda,F\right)\lesssim \frac{1}{\mu}(1-\beta){\varepsilon}\left\Vert F\right\Vert{}^{2} + \mu(1-\beta)\frac{1}{\varepsilon}{\left\Vert\varLambda\right\Vert{}_{1}^{2}}. \end{align*} Taking the parameter $$\varepsilon$$ small enough, it can be readily checked that the contributions from $$R_{j}$$ can be absorbed by Q(t) or are interpolation errors that behave as indicated in (6.7). Let us treat now the terms of the RHS of (6.8) that come from the stabilization. Using Cauchy–Schwarz’s and Young’s inequalities, we can easily control the following terms: \begin{align*} R_{9}(t) & = (1-\beta)\alpha_{u}\left(P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla{\mathscr{U}}\right],P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla E\right]\right),\\ R_{10}(t)\! & = (1-\beta)\alpha_{p}\left(P_{p}^{\perp}\left[\nabla\cdot{\mathscr{U}}\right],P_{p}^{\perp}\left[\nabla\cdot E\right]\right),\\ R_{11}(t) & \!= (1-\beta)\alpha_{u}\left(P_{u}^{\perp}\left[\nabla\mathscr{P}\right],P_{u}^{\perp}\left[\nabla G\right]\right),\\ R_{12}(t)\! & = (1-\beta)\alpha_{u}\left(P_{u}^{\perp}\left[\nabla\cdot{\varSigma}\right],P_{u}^{\perp}\left[\nabla\cdot F\right]\right),\\ R_{13}(t) & \!= \alpha_{\sigma}\!\left(P_{\sigma}^{\perp}\left[\!-(1-\!\beta)\nabla^{s}{\mathscr{U}}+\!\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla{\varSigma}\right)\right],P_{\sigma}^{\perp}\left[\!-(1-\beta)\nabla^{s}E+\!\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla F\!+g^{*}\left(\hat{u}_{h},F\right)\right)\!\right]\right)\!,\\ R_{14}(t) & \!= \alpha_{\sigma}\left(P_{\sigma}^{\perp}\left[-\frac{\lambda}{2\mu}g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right],P_{\sigma}^{\perp}\left[-(1-\beta)\nabla^{s} E+\frac{\lambda}{2\mu}\left(\hat{u}_{h}\cdot\nabla F+g^{*}\left(\hat{u}_{h},F\right)\right)\right]\right). \end{align*} Some of these terms need additional treatment. Let us consider first $$R_{9}(t)$$. When bounding it, we need to make use of the fact that \begin{align} \alpha_{u}\left\Vert P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla{\mathscr{U}}\right]\right\Vert^{2} & \leqslant \alpha_{u}\rho^{2}\left\Vert P_{u}^{\perp}\left[\left(\hat{u}_{h}-u\right)\cdot\nabla\left({\mathscr{U}-}u\right)\right]\right\Vert^{2}+\alpha_{u}\rho^{2}\left\Vert P_{u}^{\perp}\left[u\cdot\nabla\left({\mathscr{U}-}u\right)\right]\right\Vert^{2}\nonumber\\ &\quad +\alpha_{u}\rho^{2}\left\Vert P_{u}^{\perp}\left[u\cdot\nabla u\right]\right\Vert^{2}+\alpha_{u}\rho^{2}\left\Vert P_{u}^{\perp}\left[\left(\hat{u}_{h}-u\right)\cdot\nabla u\right]\right\Vert^{2}, \end{align} (6.10) where the four terms of the RHS can be shown to have structure given by (6.6) due to Assumption 6.1 and making use of adequate interpolation estimates. Terms $$R_{10}(t)$$ to $$R_{13}(t)$$ are again bounded by terms that can be absorbed by Q(t), as well as by terms that involve the orthogonal projection of an operator applied to the projection onto the finite element space of the continuous solution. They can all be treated using the triangle inequality. For example, for the term bounding $$R_{10}(t),$$ we have \begin{align} \left\Vert P_{p}^{\perp}\left[\nabla\cdot{\mathscr{U}}\right]\right\Vert^{2} \leqslant \left\Vert P_{p}^{\perp}\left[\nabla\cdot({\mathscr{U}} - u)\right]\right\Vert^{2} + \left\Vert P_{p}^{\perp}\left[\nabla\cdot u\right]\right\Vert^{2} \lesssim \Vert u \Vert^{2}_{k+1} h^{2k}, \end{align} (6.11) and can therefore be cast as (6.7). A similar strategy can be followed for the corresponding terms in $$R_{11}(t)$$, $$R_{12}(t)$$ and $$R_{13}(t)$$. The bound for this last term requires some more work. The first term can be treated as above, using the step just described for the term with $${\mathscr{U}}$$ and the same strategy as in (6.10) for $$\hat{u}_{h}\cdot \nabla{\varSigma }$$. The second term in the bound of $$R_{13}(t)$$ can be absorbed by Q(t). It only remains to treat the last term, which is particularly important because it is one of the terms responsible for the need to have $$k \geqslant d/2$$. It can be bounded as follows: \begin{align} \left\Vert P_{\sigma}^{\perp}\left[g^{*}\left(\hat{u}_{h},F\right)\right]\right\Vert^{2} &\lesssim \left\Vert P_{\sigma}^{\perp}\left[g^{*}\left(\hat{u}_{h}-u,F\right)\right]\right\Vert^{2}+\left\Vert P_{\sigma}^{\perp}\left[g^{*}\left(u,F\right)\right]\right\Vert^{2}\nonumber\\ &\lesssim \left\Vert \nabla\left(\hat{u}_{h}-u\right)\right\Vert{}^{2}\left\Vert F\right\Vert{}_{\infty}^{2}+\left\Vert \nabla u\right\Vert{}_{\infty}^{2}\left\Vert F\right\Vert{}^{2}\nonumber\\ &\lesssim \left(h^{-d}\left\Vert \nabla\left(\hat{u}_{h}-u\right)\right\Vert{}^{2}+\left\Vert \nabla u\right\Vert{}_{\infty}^{2}\right)\left\Vert F\right\Vert{}^{2}, \end{align} (6.12) where an inverse inequality has been used in the last step. Since $$\hat{U}_{h}\in \mathscr{B}_{h}$$, $$\left \Vert \nabla \left (\hat{u}_{h}-u\right )\right \Vert^{2}$$ is of order $$h^{2k}$$, and the condition $$k\geqslant d/2$$ allows us to guarantee that the term in parenthesis remains bounded as h → 0. Therefore, (6.12) (multiplied by $$\varepsilon$$) can be absorbed by Q(t). This concludes the analysis of $$R_{13}(t)$$. In the bound for $$R_{14}(t)$$, the second and third terms in the RHS have already appeared when dealing with $$R_{13}(t)$$, and thus only the first term needs to be bounded. This can be done as follows: \begin{align} \left\Vert P_{\sigma}^{\perp}\left[g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right]\right\Vert^{2} &\lesssim \left\Vert g\left(\hat{u}_{h}-u,\hat{\sigma}_{h}-\sigma\right)\right\Vert{}^{2} +\left\Vert g\left(u,\hat{\sigma}_{h}-\sigma\right)\right\Vert{}^{2} \nonumber\\ &\quad +\left\Vert P_{\sigma}^{\perp}\left[ g\left(u,\sigma\right) \right] \right\Vert^{2}+\left\Vert g\left(\hat{u}_{h}-u,\sigma\right)\right\Vert{}^{2}, \end{align} (6.13) where the four terms of the RHS can be bounded using (2.4) and inverse inequalities. Using the fact that $$\hat{U}_{h}\in \mathscr{B}_{h}$$, $$k \geqslant d/2$$ and Assumption 6.1, it follows that all terms contributing to bound $$R_{14}(t)$$ (with the adequate factor) can be written as indicated in (6.6). Term $$R_{15}(t)$$ in (6.8) corresponds to the traction or rotational term of the Oldroyd-B constitutive model, and can be written as: \begin{align} R_{15}(t) & =-\frac{\lambda}{2\mu}(g\left(u,\sigma\right),F)+\frac{\lambda}{2\mu}(g\left(\hat{u}_{h},\hat{\sigma}_{h}\right),F)\nonumber\\ & = -\frac{\lambda}{2\mu}(g\left(u,\sigma-\hat{\sigma}_{h}\right),F)+\frac{\lambda}{2\mu}(g\left(u-\hat{u}_{h},\sigma-\hat{\sigma}_{h}\right),F) -\frac{\lambda}{2\mu}(g\left(u-\hat{u}_{h},\sigma\right),F). \end{align} (6.14) The three terms in the RHS can be bounded using Cauchy–Schwarz’s and Young’s inequalities and the inverse inequality (2.2). When these inequalities are used in (6.14) it is seen that we again recover (6.6). The remaining $$R_{i}$$ terms are all multiplied by $$\gamma$$. To distinguish these terms from the others we use $$\xi$$ instead of $$\varepsilon$$. The first of them can be bounded using Cauchy–Schwarz’s and Young’s inequalities, as well as the inverse inequality, as $$R_{16}(t) \leqslant (1-\beta)\gamma\alpha_{u}\left\Vert \nabla\varPi\right\Vert \left\Vert{P_{u,0}\left[\nabla G\right]}\right\Vert \lesssim \frac{1}{\mu}(1-\beta)\gamma\frac{1}{\xi}\left\Vert \varPi\right\Vert{}^{2} + (1-\beta)\gamma\mu h^{-2} {\alpha_{u}^{2}} \xi\left\Vert{P_{u,0}\left[\nabla G\right]}\right\Vert{}^{2}.$$ The next term requires some more elaboration. Observe first that $$( \partial_{t}(u-u_{h}),{P_{u,0}\left[\nabla G\right]}) = ( \partial_{t} u,{P_{u,0}\left[\nabla G\right]}) + ( \partial_{t} \nabla\cdot u_{h}, G).$$ Since u is divergence free and vanishes on $$\partial \varOmega$$, we have that $$( \partial_{t} u,{P_{u,0}\left[\nabla G\right]}) = ( \partial_{t} u, P_{u}[\nabla G] - \nabla G) = -\left(P_{u}^{\bot}[\partial_{t} u], P_{u}^{\bot}[\nabla G]\right).$$ Making use of the continuity equation (5.2), we can write $$( \partial_{t} \nabla\cdot u_{h}, G) = -\,\alpha_{u} \left(\partial_{t} P_{u}^{\bot} [\nabla p_{h}], P_{u}^{\bot} [\nabla G]\right) = \alpha_{u} \left(\partial_{t} P_{u}^{\bot} [\nabla G], P_{u}^{\bot} [\nabla G]\right) - \alpha_{u} \left(\partial_{t} P_{u}^{\bot} [\nabla{\mathscr{P}} ], P_{u}^{\bot} [\nabla G]\right).$$ From these two equalities, we obtain \begin{align*} R_{17}(t) =& \,(1-\beta)\gamma\alpha_{u}\rho ( \partial_{t}(u-u_{h}),{P_{u,0}\left[\nabla G\right]}) \\ =& -(1-\beta)\gamma\alpha_{u}\rho \left( P_{u}^{\bot}[\partial_{t} u], P_{u}^{\bot} [\nabla G]\right)\\ & + (1-\beta){\gamma\alpha_{u}^{2}} \rho \frac{1}{2}\frac{\mathrm{d}}{\mathrm{d} t} \left\Vert P_{u}^{\bot} [\nabla G] \right\Vert^{2} - (1-\beta){\gamma\alpha_{u}^{2}} \rho \left(\partial_{t} P_{u}^{\bot} [\nabla{\mathscr{P}}], P_{u}^{\bot} [\nabla G]\right). \end{align*} After using Schwarz’s and Young’s inequalities, the first term in the RHS can be treated as $$R_{4}(t)$$ (note that $$\varLambda = P_{u}^{\bot } [u]$$), and the third term can be treated as $$R_{11}(t)$$; therefore, they both can be cast as (6.7). It only remains to deal with the second term. Let $$\chi_{h} (G(t)) := (1-\beta)\alpha_{u} \left\Vert P_{u}^{\bot} [\nabla G(t)] \right\Vert^{2} \geqslant 0.$$ After integrating in time, in (6.9) we have the terms \begin{align} {\int_{0}^{T}} \chi_{h} (G(t))\,\mathrm{d} t + \dots \lesssim \alpha_{u} \rho \gamma \chi_{h}(G(T)) + \dots, \end{align} (6.15) the one on the RHS coming from $$R_{17}(t)$$. Since $$\hat{U}_{h}\in \mathscr{B}_{h}$$ and because of the boundedness of U described in Assumption 6.1, the RHS in (5.14) is bounded for all t. Since, as a consequence of the proof of Theorem 5.2, $$p_{h}$$ is continuous in $$[0,\infty )$$, we have that $$\alpha _{u} \Vert P_{u}^{\bot } [\nabla p_{h}(t)]\Vert ^{2}$$ is continuous and bounded in t. Assumption 6.1 implies that $$\alpha _{u} \Vert P_{u}^{\bot } [\nabla p(t)]\Vert ^{2}$$ is also continuous and bounded in t, and therefore $$\chi _{h} (G(T))$$ is continuous and bounded for all T. In this case, there exists $$T_{0}$$ such that \begin{align} \alpha_{u} \rho \chi_{h}(G(T)) \leqslant{\int_{0}^{T}} \chi_{h}(G(t))\,\mathrm{d}t =:\mathscr{G}(T), \quad \textrm{for all}~T \geqslant T_{0}. \end{align} (6.16) Let us prove this. The LHS is bounded in T, and thus the result is obvious if $$\mathscr{G}(T)$$ is unbounded. Suppose now that $$\mathscr{G}(T) \nearrow M < \infty$$ as $$T \to \infty$$. In this case, $$\chi _{h}(G(t)) \to 0^{+}$$ as $$t\ \to\ \infty$$. Thefore, there must exist $$T_{1}$$ such that $$\mathscr{G}(T)> M/2$$ for all $$T> T_{1}$$ and $$T_{0}> T_{1}$$ such that $$\alpha _{u} \rho \chi _{h}(G(T)) < M/2$$ for all $$T> T_{0}$$, so that (6.16) holds. Once (6.16) is proved, it is observed that the term in the RHS of (6.15) can be absorbed by the one in the LHS for an appropriate $$\gamma$$. This concludes the analysis of $$R_{17}(t)$$. The following terms can be controlled using the inverse inequality and Young’s inequality: \begin{align*} R_{18}(t) & \lesssim \beta(1-\beta)\gamma\alpha_{u}\mu\left(\left\Vert \nabla^{s}E\right\Vert +\left\Vert \nabla^{s}\varLambda\right\Vert \right)\left\Vert \nabla^{s}\left({P_{u,0}\left[\nabla G\right]}\right)\right\Vert\!,\\ R_{19}(t) & \leqslant (1-\beta)\gamma\alpha_{u}\left(\left\Vert F\right\Vert +\left\Vert \varGamma\right\Vert \right)\left\Vert \nabla^{s}\left({P_{u,0}\left[\nabla G\right]}\right)\right\Vert\!. \end{align*} To see that the bounds obtained fall within the structure (6.6), we have to use the expression of the stabilization parameter $$\alpha _{u}$$ to check that $$\mu h^{-2}\lesssim \alpha _{u}^{-1}$$ and take $$\gamma$$ small enough to balance $$\xi ^{-1}$$, which may need to be large. The next terms to consider come from the stabilization \begin{align*} & R_{20}(t) =\!(1-\beta){\gamma\alpha_{u}^{2}}\left(P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla\left(u_{h}-{\mathscr{U}} +{\mathscr{U}}\right)\right],P_{u}^{\perp}\left[\rho\hat{u}_{h}\cdot\nabla\left({P_{u,0}\left[\nabla G\right]}\right)\right]\right)\!,\\ & R_{21}(t) =\!(1-\beta)\gamma\alpha_{u}\alpha_{p}\left(P_{p}^{\perp}\left[\nabla\cdot\left(u_{h}-{\mathscr{U}} +{\mathscr{U}}\right)\right],P_{p}^{\perp}\left[\nabla\cdot\left({P_{u,0}\left[\nabla G\right]}\right)\right]\right)\!,\\ & R_{22}(t) =\!(1\!-\!\beta)\gamma\alpha_{u}\alpha_{\sigma}\!\left(\!P_{\sigma}^{\perp}\left[\!(1\!-\!\beta)\nabla^{s}\!\left(u_{h}\!-{\mathscr{U}} \!+{\mathscr{U}}\right)\!-\!\frac{\lambda}{2\mu}\hat{u}_{h}\cdot\nabla\left(\sigma_{h}\!-\!{\varSigma} \!+\!{\varSigma}\right)\!\right]\!, P_{\sigma}^{\perp}\left[\nabla^{s}\!\left({P_{u,0}\left[\nabla G\right]}\right)\right]\right)\!,\\ & R_{23}(t) =\!(1-\beta)\gamma\alpha_{u}\alpha_{\sigma}\left(P_{\sigma}^{\perp}\left[\frac{\lambda}{2\mu}g\left(\hat{u}_{h},\hat{\sigma}_{h}\right)\right],P_{\sigma}^{\perp}\left[\nabla^{s}\left({P_{u,0}\left[\nabla G\right]}\right)\right]\right). \end{align*} All these terms can be bounded using a similar strategy as before. With this we conclude the analysis of the last term in (6.8) (or in (6.9)). The second term, S(t), is easy to treat; after using Young’s inequality, the first term that arises can be absorbed by Q(t) and the second bounded as in (6.12). We have written S(t) independently to emphasize that it arises because of the linearization of the rotational term: since we evaluate it with $$\hat{u}_{h}$$ and $$\hat{\sigma }_{h}$$, it does not contribute to stability (that is to say, its contribution does not appear in Q(t)). We could however have considered it as one more $$R_{j}(t)$$ term. It only remains to deal with the convective terms N(t) in (6.8). The first two of these terms can be written as follows: \begin{align} & \rho(1-\beta)\left(\tilde{c}\left(u,u,E\right)-\tilde{c}\left(\hat{u}_{h},u_{h},E\right)\right)\nonumber \\ & \quad = \rho(1-\beta)\left(\tilde{c}\left(\left(u-\hat{u}_{h}\right),u,E\right)-\tilde{c}\left(\left(u-\hat{u}_{h}\right),\varLambda,E\right)+\tilde{c}\left(u,\varLambda,E\right)\right)\!. \end{align} (6.17) The three RHS terms of the above equation can be easily bounded fitting the structure (6.6)–(6.7). A similar procedure can be applied to the convective terms arising from the constitutive equation, which can be written as: $$\frac{\lambda}{2\mu}\left(\tilde{c}\left(u,\sigma,F\right)-\tilde{c}\left(\hat{u}_{h},\sigma_{h},F\right)\right) = \frac{\lambda}{2\mu}\left(\tilde{c}\left(\left(u-\hat{u}_{h}\right),\sigma,F\right)-\tilde{c}\left(\left(u-\hat{u}_{h}\right),\varGamma,F\right)+\tilde{c}\left(u,\varGamma,F\right)\right)\!.$$ The remaining convective terms to be controlled in the LHS of (6.8) can be written as \begin{align*} &(1-\beta)\gamma\alpha_{u}\rho\left(\tilde{c}\left(u,u,{P_{u,0}\left[\nabla G\right]}\right)-\tilde{c}\left(\hat{u}_{h},u_{h},{P_{u,0}\left[\nabla G\right]}\right)\right)\\ & \quad = (1-\beta)\gamma\alpha_{u}\rho\left( \tilde{c}(u-\hat{u}_{h},u,{P_{u,0}\left[\nabla G\right]}) - \tilde{c}(u-\hat{u}_{h},u -u_{h},{P_{u,0}\left[\nabla G\right]}) + \tilde{c}(u,u-u_{h},{P_{u,0}\left[\nabla G\right]}) \right)\!. \end{align*} The strategy to bound these terms is similar to what we have used heretofore. This concludes the proof of the third step. Step 4. Fixed point theorem. According to Step 3, the ball $$\mathscr{B}_{h}$$ is invariant under the map $$\delta$$, that is to say, $$\delta \left (\mathscr{B}_{h}\right )\subset \mathscr{B}_{h}$$. Therefore, applying Schauder’s fixed point theorem, we conclude that there exists $$(u_{h},p_{h},\sigma _{h})\in \mathscr{B}_{h}$$ such that $$(u_{h},p_{h},\sigma _{h}) = \delta (u_{h},p_{h},\sigma _{h})$$, which, in view of the definition of $$\delta$$, implies that $$(u_{h},p_{h},\sigma _{h})$$ is a solution of (4.1). Because of the definition of $$\mathscr{B}_{h}$$ in (6.2), we also obtain an error estimate. 7. Conclusions In this article we have presented the numerical analysis of a stabilized finite element approximation proposed to solve viscoelastic fluid flows, in the nonlinear time-dependent case. This analysis has confirmed the numerical results obtained in other works, where the method was proposed and tested in nonlinear examples. In particular, we have shown this using a fixed point theorem under suitable regularity conditions in velocity, stress and pressure. As it is usual for nonlinear problems involving the Navier–Stokes equations, the estimates do not show information about their behaviour when the local Reynolds and Weissenberg numbers increase. On the other hand, some control over the pressure, although very weak, has been obtained. Funding E.C. thanks the support given by the Chilean Council for Scientific and Technological Research (CONICYT-FONDECYT 11160160), and from the Universidad de S. de Chile. R.C. gratefully acknowledges the support received from the Catalan Government through the Catalan Institution for Research and Advanced Studies (ICREA) Acadèmia Research Program. References Baaijens , F . ( 1998 ) Mixed finite element methods for viscoelastic flow analysis: a review . J. Nonnewton. Fluid Mech. , 79 , 361 – 385 . Google Scholar CrossRef Search ADS Badia , S. & Codina , R. ( 2009 ) On a multiscale approach to the transient Stokes problem: dynamic subscales and anisotropic space–time discretization . Appl. Math. Comput. , 207 , 415 – 433 . Badia , S. , Codina , R. & Gutiérrez-Santacreu , J. V. ( 2010 ) Long-term stability estimates and existence of a global attractor in a finite element approximation of the Navier–Stokes equations with numerical subgrid scale modeling . SIAM J. Numer. Anal. , 48 , 1013 – 1037 . Google Scholar CrossRef Search ADS Baranger , J. & Sandri , D. ( 1992 ) Finite element approximation of viscoelastic fluid flow: existence of approximate solutions and error bounds . Numer. Math. , 63 , 13 – 27 . Google Scholar CrossRef Search ADS Baranger , J. & Wardi , S. ( 1995 ) Numerical analysis of a FEM for a transient viscoelastic flow . Comput. Methods Appl. Mech. Eng. , 125 , 171 – 185 . Google Scholar CrossRef Search ADS Barrett , J. & Boyaval , S. ( 2011 ) Existence and approximation of a (regularized) oldroyd-b model . Math. Models Meth. Appl. Sci. , 21 , 1783 – 1837 . Google Scholar CrossRef Search ADS Bird , R. B. , Armstrong , R. C. & Hassager , O. ( 1987a ) Dynamics of Polymeric Liquids, vol. 1: Fluid Mechanics , 2nd edn. New York : Wiley . Bird , R. B. , Armstrong , R. C. & Hassager , O. ( 1987b ) Dynamics of Polymeric Liquids, vol. 2: Kinetic Theory , 2nd edn. New York : Wiley . Bird , R. B. , Stewart , W. E. & Lightfoot , E. N. ( 2002 ) Transport Phenomena . New York : Wiley . Bogaerds , A. , Verbeeten , W. , Peters , G. & Baaijens , F. ( 1999 ) 3D viscoelastic analysis of a polymer solution in a complex flow . Comput. Meth. Appl. Mech. Eng. , 180 , 413 – 430 . Google Scholar CrossRef Search ADS Bonito , A. & Burman , E. ( 2007 ) A continuous interior penalty method for viscoelastic flows . SIAM J. Sci. Comput. , 30 , 1156 – 1177 . Google Scholar CrossRef Search ADS Bonito , A. , Clément , P. & Picasso , M. ( 2007 ) Mathematical and numerical analysis of a simplified time-dependent viscoelastic flow . Num. Math. , 107 , 213 – 255 . Google Scholar CrossRef Search ADS Burman , E. & Fernández , M. ( 2007 ) Continuous interior penalty finite element method for the time-dependent Navier–Stokes equations: space discretization and convergence . Num. Math. , 107 , 39 – 77 . Google Scholar CrossRef Search ADS Castillo , E. , Baiges , J. & Codina , R. ( 2015 ) Approximation of the two-fluid flow problem for viscoelastic fluids using the level set method and pressure enriched finite element shape functions . J. Nonnewton. Fluid Mech. , 225 , 37 – 53 . Google Scholar CrossRef Search ADS Castillo , E. & Codina , R. ( 2014a ) Stabilized stress-velocity-pressure finite element formulations of the Navier-Stokes problem for fluids with non-linear viscosity . Comput. Methods Appl. Mech. Eng. , 279 , 554 – 578 . Google Scholar CrossRef Search ADS Castillo , E. & Codina , R. ( 2014b ) Variational multi-scale stabilized formulations for the stationary three-field incompressible viscoelastic flow problem . Comput. Methods Appl. Mech. Eng. , 279 , 579 – 605 . Google Scholar CrossRef Search ADS Castillo , E. & Codina , R. ( 2015 ) First, second and third order fractional step methods for the three-field viscoelastic flow problem . J. Comput. Phys. , 296 , 113 – 137 . Google Scholar CrossRef Search ADS Castillo , E. & Codina , R. ( 2017a ) Finite element approximation of the viscoelastic flow problem: a non-residual based stabilized formulation . Comput. Fluids , 142 , 72 – 78 . Google Scholar CrossRef Search ADS Castillo , E. & Codina , R. ( 2017b ) Numerical analysis of a stabilized finite element approximation for the three-field linearized viscoelastic fluid problem using arbitrary interpolations . ESAIM Math. Model. Num. Anal. , 51 , 1407 – 1427 . Ciarlet , P. ( 2013 ) Linear and Nonlinear Functional Analysis with Applications . Philadelphia, PA, USA : Society for Industrial and Applied Mathematics . Codina , R . ( 2000 ) Stabilization of incompressibility and convection through orthogonal sub-scales in finite element methods . Comput. Methods Appl. Mech. Eng. , 190 , 1579 – 1599 . Google Scholar CrossRef Search ADS Codina , R . ( 2009 ) Finite element approximation of the three-field formulation of the Stokes problem using arbitrary interpolations . SIAM J. Numer. Anal. , 47 , 699 – 718 . Google Scholar CrossRef Search ADS Codina , R. & Blasco , J. ( 2000 ) Stabilized finite element method for the transient Navier–Stokes equations based on a pressure gradient projection . Comput. Methods Appl. Mech. Eng. , 182 , 277 – 300 . Google Scholar CrossRef Search ADS Codina , R. , Principe , J. , Guasch , O. & Badia , S. ( 2007 ) Time dependent subscales in the stabilized finite element approximation of incompressible flow problems . Comput. Methods Appl. Mech. Eng. , 196 , 2413 – 2430 . Google Scholar CrossRef Search ADS Ern , A. & Guermond , J.-L. ( 2004 ) Theory and Practice of Finite Elements , vol. 159, 1st edn. New York : Springer . Google Scholar CrossRef Search ADS Ervin , V. , Lee , H. & Ntasin , L. ( 2005 ) Analysis of the Oseen-viscoelastic fluid flow problem . J. Nonnewton. Fluid Mech. , 127 , 157 – 168 . Google Scholar CrossRef Search ADS Ervin , V. & Miles , W. ( 2003 ) Approximation of time-dependent viscoelastic fluid flow: SUPG approximation . SIAM J. Numer. Anal. , 41 , 457 – 486 . Google Scholar CrossRef Search ADS Ervin , V. & Miles , W. ( 2005 ) Approximation of time-dependent, multi-component, viscoelastic fluid flow . Comput. Methods Appl. Mech. Eng. , 194 , 2229 – 2255 . Google Scholar CrossRef Search ADS Fernández-Cara , E. , Guillén , F. & Ortega , R. ( 2002 ) Mathematical modeling and analysis of viscoelastic fluids of the Oldroyd kind . Handbook of Numerical Analysis , vol. 8. Amsterdam : North-Holland , pp. 543 – 661 . Fortin , M. & Fortin , A. ( 1989 ) A new approach for the FEM simulation of viscoelastic flows . J. Nonnewton. Fluid Mech. , 32 , 295 – 310 . Google Scholar CrossRef Search ADS Girault , V. & Raviart , P.-A. ( 1986 ) Finite Elements for the Navier–Stokes Equations . Berlin : Springer . Google Scholar CrossRef Search ADS Guillopé , C. & Saut , J. ( 1990 ) Existence results for the flow of viscoelastic fluids with a differential constitutive law . Nonlinear Anal. Theory Methods Appl. , 15 , 849 – 869 . Google Scholar CrossRef Search ADS Howell , J . ( 2009 ) Dual-mixed finite element approximation of Stokes and nonlinear Stokes problems using trace-free velocity gradients . J. Comput. Appl. Math. , 231 , 780 – 792 . Google Scholar CrossRef Search ADS Hughes , T. , Feijóo , G. , Mazzei , L. & Quincy , J.-B. ( 1998 ) The variational multiscale method-a paradigm for computational mechanics . Comput. Methods Appl. Mech. Eng. , 166 , 3 – 24 . Google Scholar CrossRef Search ADS Lions , P. L. & Masmoudi , N. ( 2000 ) Global solutions for some Oldroyd models of non-Newtonian flows . Chin. Ann. Math. Ser. B , 21 , 131 – 146 . Google Scholar CrossRef Search ADS Marchal , J. & Crochet , M. ( 1987 ) A new mixed finite element for calculating viscoelastic flow . J. Nonnewton. Fluid Mech. , 26 , 77 – 114 . Google Scholar CrossRef Search ADS Pani , A. K. & Yuan , J. Y. ( 2005 ) Semidiscrete finite element Galerkin approximations to the equations of motion arising in the Oldroyd model . IMA J. Numer. Anal. , 25 , 750 – 782 . Google Scholar CrossRef Search ADS Picasso , M. & Rappaz , J. ( 2001 ) Existence, a priori and a posteriori error estimates for a nonlinear three-field problem arising from Oldroyd-B viscoelastic flows . ESAIM: Math. Model. Numer. Anal. , 35 , 879 – 897 . Google Scholar CrossRef Search ADS Renardy , M . ( 1985 ) Existence of slow steady flows of viscoelastic fluids with differential constitutive equations . ZAMM/J. Appl. Math. Mech. , 65 , 449 – 451 . Google Scholar CrossRef Search ADS Renardy , M. ( 2000 ) Mathematical analysis of viscoelastic flows . CBMS-NSF Regional Conference Series in Applied Mathematics . SIAM : Philadelphia . Google Scholar CrossRef Search ADS Sandri , D . ( 1994 ) Finite element approximation of viscoelastic fluid flow: existence of approximate solutions and error bounds. Continuous approximation of the stress . SIAM J. Numer. Anal. , 31 , 362 – 377 . Google Scholar CrossRef Search ADS Sureshkumar , R. & Beris , A. ( 1995 ) Effect of artificial stress diffusivity on the stability of numerical calculations and the flow dynamics of time-dependent viscoelastic flows . J. Nonnewton. Fluid Mech. , 60 , 53 – 80 . Google Scholar CrossRef Search ADS © The Author(s) 2018. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) For permissions, please e-mail: journals. permissions@oup.com

### Journal

IMA Journal of Numerical AnalysisOxford University Press

Published: Apr 24, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations