John W Barrett

S Ébastien Boyaval

FINITE ELEMENT APPROXIMATION OF THE FENE-P MODEL

Keywords: FENE-P model, entropy, finite element method, convergence analysis, stress diffusion, existence of weak solutions AMS Subject Classification: 35Q30, 65M12, 65M60, 76A10, 76M10, 82D60

We extend our analysis on the Oldroyd-B model in Barrett and Boyaval [1] to consider the finite element approximation of the FENE-P system of equations, which models a dilute polymeric fluid, in a bounded domain D ⊂ R d , d = 2 or 3, subject to no flow boundary conditions. Our schemes are based on approximating the pressure and the symmetric conformation tensor by either (a) piecewise constants or (b) continuous piecewise linears. In case (a) the velocity field is approximated by continuous piecewise quadratics (d = 2) or a reduced version, where the tangential component on each simplicial edge (d = 2) or face (d = 3) is linear. In case (b) the velocity field is approximated by continuous piecewise quadratics or the mini-element. We show that both of these types of schemes, based on the backward Euler type time discretization, satisfy a free energy bound, which involves the logarithm of both the conformation tensor and a linear function of its trace, without any constraint on the time step. Furthermore, for our approximation (b) in the presence of an additional dissipative term in the stress equation, the so-called FENE-P model with stress diffusion, we show (subsequence) convergence in the case d = 2, as the spatial and temporal discretization parameters tend to zero, towards global-intime weak solutions of this FENE-P system. Hence, we prove existence of global-in-time weak solutions to the FENE-P model with stress diffusion in two spatial dimensions.

Here R d×d S denotes the set of symmetric R d×d matrices, and R d×d S,>0 the set of symmetric positive definite R d×d matrices. In addition, I ∈ R d×d S,>0 is the identity, and tr(•) denotes trace. The unknowns in (P) are the velocity of the fluid, u, the hydrostatic pressure, p, and the symmetric conformation tensor of the polymer molecules, σ. The latter is linked to the symmetric polymeric extra-stress tensor τ through the relation τ = ε Wi A(σ) σ. In addition, f : (t, x) ∈ D T → f (t, x) ∈ R d is the given density of body forces acting on the fluid; and the following given parameters are dimensionless: the Reynolds number Re ∈ R >0 , the Weissenberg number Wi ∈ R >0 , the elasticto-viscous viscosity fraction ε ∈ (0, 1), and the FENE-P parameter b > 0 (related to a maximal admissible extensibility of the polymer molecules within the fluid). For the sake of simplicity, we will limit ourselves to the no flow boundary conditions (1.1f). Finally, we denote ∇u(t, x) ∈ R d×d the velocity gradient tensor field with [∇u] ij = ∂ui ∂xj , and (div σ)(t, x) ∈ R d the vector field with [div σ] i = d j=1 ∂σij ∂xj . For data f ≡ 0, divergence free u 0 ∈ [L 2 (D)] d , and σ 0 , which is symmetric positive definite for a.e. x ∈ D, satisfying ln(1 -tr(σ 0 ) b ) ∈ L 1 (D), then the existence of a global-in-time weak solution u ∈ L ∞ (0, T ; [L 2

(D)] d ) ∩ L 2 (0, T, [H 1 0 (D)] d ), σ ∈ L ∞ (0, T ; [L ∞ (D)] d×d
) and τ ∈ L 2 (0, T ; [L 2 (D)] d×d ) to (P), (1.1a-f), was proved in Masmoudi [START_REF] Masmoudi | Global existence of weak solutions to macroscopic models of polymeric flows[END_REF].

In this work, we consider finite element approximations of the FENE-P system (P) and the corresponding model with stress diffusion, (P α ), which is obtained by adding the dissipative term α ∆σ for a given α ∈ R >0 to the right-hand side of (1.1c) with an additional no flux boundary condition for σ on ∂D. This paper extends the results in Barrett and Boyaval [START_REF] Barrett | Existence and approximation of a (regularized) Oldroyd-B model[END_REF], where finite element approximations of the corresponding Oldroyd-B models, where A(σ) = Iσ -1 , were introduced and analysed. In fact, the convergence proof of the finite element approximation of the Oldroyd-B model with stress diffusion for d = 2 in [START_REF] Barrett | Existence and approximation of a (regularized) Oldroyd-B model[END_REF] provided the first existence proof of global-in-time weak solutions for this system. Note that A(σ) = Iσ -1 is the formal limit of (1.2) for infinite extensibility; that is, b → ∞.

The model (P α ) has been considered computationally in Sureshkumar and Beris [START_REF] Sureshkumar | Effect of artificial stress diffusivity on the stability of numerical calculations and the flow dynamics of time-dependent viscoelastic flows[END_REF]. We recall also that El-Kareh and Leal [START_REF] El-Kareh | Existence of solutions for all Deborah numbers for a non-Newtonian model modified to include diffusion[END_REF] showed the existence of a weak solution to a modified stationary FENE-P system of equations, which included stress diffusion, but there an additional regularization was also present in their modified system and played an essential role in their proof. We stress that the dissipative term α ∆σ in (P α ) is not a regularization, but can be physically motivated through the centre-of-mass diffusion in the related microscopic-macroscopic polymer model, though with a positive α ≪ 1, see Barrett and Süli [START_REF] Barrett | Existence of global weak solutions to some regularized kinetic models of dilute polymers[END_REF], [START_REF] Barrett | Existence of global weak solutions to finitely extensible nonlinear bead-spring chain models for dilute polymers with variable density and viscosity[END_REF], Schieber [START_REF] Schieber | Generalized Brownian configuration field for Fokker-Planck equations including center-of-mass diffusion[END_REF] and Degond and Liu [START_REF] Degond | Kinetic models for polymers with inertial effects[END_REF].

Barrett and Süli have introduced, and proved the existence of global-in-time weak solutions for d = 2 and 3 to, microscopic-macroscopic dumbbell models of dilute polymers with center-of-mass diffusion in the corresponding Fokker-Planck equation for a finitely extensible nonlinear elastic (FENE) spring law or a Hookean-type spring law, see [START_REF] Barrett | Existence and equilibration of global weak solutions to kinetic models for dilute polymers I: Finitely extensible nonlinear bead-spring chains[END_REF] and [START_REF] Barrett | Existence and equilibration of global weak solutions to kinetic models for dilute polymers II: Hookean-type bead-spring chains[END_REF]. Recently, Barrett and Süli [START_REF] Barrett | Existence of global weak solutions to the kinetic Hookean dumbbell model for incompressible dilute polymeric fluids[END_REF] have proved rigorously that the macroscopic Oldroyd-B model with stress diffusion is the exact closure of the microscopic-macroscopic Hookean dumbbell model with center-of-mass diffusion for d = 2, when the existence of global-in-time weak solutions to both models can be proved. In addition, Barrett and Süli [START_REF] Barrett | Finite element approximation of finitely extensible nonlinear elastic dumbbell models for dilute polymers[END_REF] have introduced and analysed a finite element approximation for the FENE microscopic-macroscopic dumbbell model with center-of-mass diffusion.

From a physical viewpoint, the FENE-P model is more realistic than the Oldroyd-B model because it accounts for the finite-extensibility of the polymer molecules in the fluid through the non-dimensional parameter b > 0. From a mathematical viewpoint, compared to the Oldroyd-B model where the nonlinear terms are only the material derivative terms (like (∇u)σ), the FENE-P model has an additional singular nonlinearity due to the factor 1 -tr(•) b -1 in the definition of A(•), which necessitates a careful mathematical treatment. Hence, this paper is not a trivial extension of [START_REF] Barrett | Existence and approximation of a (regularized) Oldroyd-B model[END_REF]. In fact, the latter additional nonlinearity is exactly what makes the FENE-P model closer to the physics of polymers than the Oldroyd-B model, and thus also to many other macroscopic models based on different constitutive relations that have been developed by physicists for polymers. We note the FENE-P system is the approximate macroscopic closure of the FENE microscopic-macroscopic dumbbell model, whereas the Oldroyd-B system is the exact macroscopic closure of the Hookean microscopic-macroscopic dumbbell model. Hence, the microscopic-macroscopic dumbbell models corresponding to Oldroyd-B and FENE-P, only the spring laws differ; see e.g. Bird et al. [START_REF] Bird | Dynamics of Polymeric Liquids[END_REF] and Renardy [START_REF] Renardy | Mathematical Analysis of Viscoelastic Flows[END_REF] for a more complete review of the differences between the Oldroyd-B and the FENE-P models from the physical viewpoint, and for other macroscopic models with more nonlinear effects than the Oldroyd-B model, e.g. the Giesekus model and the Phan-Thien Tanner model.

Similarly to [START_REF] Barrett | Existence and approximation of a (regularized) Oldroyd-B model[END_REF], our analysis in the present paper exploits the underlying free energy of the system, see Wapperom and Hulsen [START_REF] Wapperom | Thermodynamics of viscoelastic fluids: The temperature equation[END_REF] and Hu and Lelièvre [START_REF] Hu | New entropy estimates for the Oldroyd-B model, and related models[END_REF]. In particular, the finite element approximation of (P α ) has to be constructed extremely carefully to inherit this free energy structure, and requires the approximation of tr(σ) as a new unknown. It is definitely not our goal to review all the macroscopic models used in rheology, although similar studies could probably be pursued for other macroscopic models endowed with a free energy. We will point out the main differences with Barrett and Boyaval [START_REF] Barrett | Existence and approximation of a (regularized) Oldroyd-B model[END_REF], and we thus hope to sufficiently suggest how our technique could be adapted to any nonlinear model with a free energy. We believe that our approach contributes to a better understanding of the numerical stability of the models used in computational rheology, where numerical instabilities sometimes termed "High-Weissenberg Number Problems", see HWNP in Owens and Phillips [START_REF] Owens | Computational Rheology[END_REF], still persist. Indeed, as exposed in Boyaval et al. [START_REF] Boyaval | Free-energy-dissipative schemes for the Oldroyd-B model[END_REF], our point is that to make progress in this area one should identify sufficiently general rules for the derivation of good discretizations of macroscopic models such that they retain the dissipative structure of weak solutions to the system, at least in some benchmark flows.

The outline of this paper is as follows. First, we end this section by introducing our notation and some auxiliary results. In Section 2 we review the formal free energy bound for the FENE-P system (P). In Section 3 we introduce our regularization G δ of of G ≡ ln, which appears in the definition of the free energy of the FENE-P system (P). We then introduce a regularized problem (P δ ), and show a formal free energy bound for it. In Section 4, on assuming that D is a polytope for ease of exposition, we introduce our finite element approximation of (P δ ), namely (P ∆t δ,h ), based on approximating the pressure and the symmetric conformation tensor by piecewise constants; and the velocity field with continuous piecewise quadratics or a reduced version, where the tangential component on each simplicial edge (d = 2) or face (d = 3) is linear. Using the Brouwer fixed point theorem, we prove existence of a solution to (P ∆t δ,h ) and show that it satisfies a discrete regularized free energy bound for any choice of time step; see Theorem 4.1. We conclude by showing that, in the limit δ → 0 + , these solutions of (P ∆t δ,h ) converge to a solution of (P ∆t h ) with the approximation of the conformation tensor being positive definite and having a trace strictly less than b. Moreover, this solution of (P ∆t h ) satisfies a discrete free energy bound; see Theorem 4.2. Next, in Section 5 we introduce the FENE-P system with stress diffusion, (P α ), where the dissipative term α ∆σ has been added to the right-hand side of (1.1c). We then introduce the corresponding regularized version (P α,δ ), and show a formal free energy bound for it. In Section 6 we introduce our finite element approximation of (P α,δ ), namely (P ∆t α,δ,h ), based on approximating the velocity field with continuous piecewise quadratics or the mini element, and the pressure, the symmetric conformation tensor and its trace by continuous piecewise linears. Here we assume that D is a convex polytope and that the finite element mesh consists of quasi-uniform non-obtuse simplices. Using the Brouwer fixed point theorem, we prove existence of a solution to (P ∆t α,δ,h ) and show that it satisfies a discrete regularized free energy bound for any choice of time step; see Theorem 6.1. In Section 7 we prove, in the case d = 2, (subsequence) convergence of the solutions of (P ∆t α,δ,h ), as the regularization parameter, δ, and the spatial, h, and temporal, ∆t, discretization parameters tend to zero, to global-in-time weak solutions of (P α ); see Theorem 7.3. This existence result for (P α ) is new to the literature. 1.2. Notation and auxiliary results. The absolute value and the negative part of a real number s ∈ R are denoted by |s| := max{s, -s} and [s] -= min{s, 0}, respectively. We adopt the following notation for inner products

v • w := d i=1 v i w i ≡ v T w = w T v ∀v, w ∈ R d , (1.3a) φ : ψ := d i=1 d j=1 φ ij ψ ij ≡ tr φ T ψ = tr ψ T φ ∀φ, ψ ∈ R d×d , (1.3b) ∇φ :: ∇ψ := d i=1 d j=1 ∇φ ij • ∇ψ ij ∀φ, ψ ∈ R d×d ; (1.3c)
where • T and tr (•) denote transposition and trace, respectively. The corresponding norms are

v := (v • v) 1 2 , ∇v := (∇v : ∇v) 1 2 ∀v ∈ R d ; (1.4a) φ := (φ : φ) 1 2 ,
∇φ := (∇φ :: ∇φ)

1 2 ∀φ ∈ R d×d . (1.4b)
We will use on several occasions that tr(φ) = tr(φ T ) and tr(φψ) = tr(ψφ) for all φ, ψ ∈ R d×d , and

φχ T : ψ = χφ : ψ = χ : ψφ ∀φ, ψ ∈ R d×d S , χ ∈ R d×d , (1.5a) ψφ ≤ ψ φ ∀φ, ψ ∈ R d×d , (1.5b) φ v ≤ φ v ∀φ ∈ R d×d , v ∈ R d . (1.5c) For any φ ∈ R d×d S
, there exists a decomposition

(1.6) φ = O T DO ⇒ tr (φ) = tr (D) ,
where O ∈ R d×d is an orthogonal matrix and D ∈ R d×d a diagonal matrix. Hence, for any

g : R → R, one can define g(φ) ∈ R d×d S as (1.7) g(φ) := O T g(D)O ⇒ tr (g(φ)) = tr (g(D)) , where g(D) ∈ R d×d S is the diagonal matrix with entries [g(D)] ii = g(D ii ), i = 1, . . . , d.
Although the diagonal decomposition (1.6) is not unique, (1.7) uniquely defines g(φ). We note for later purposes that

d -1 (tr(|φ|)) 2 ≤ φ 2 ≤ (tr(|φ|)) 2 ∀φ ∈ R d×d S . (1.8)
One can show via diagonalization, see e.g. [START_REF] Barrett | Existence and approximation of a (regularized) Oldroyd-B model[END_REF] for details, that for all concave function g ∈ C 1 (R), it holds

(1.9) (φ -ψ) : g ′ (ψ) ≥ tr (g(φ) -g(ψ)) ≥ (φ -ψ) : g ′ (φ) ∀φ, ψ ∈ R d×d S ,
where g ′ denotes the first derivative of g. If g ∈ C 1 (R) is convex, the inequalities in (1.9) are reversed. It follows from (1.9) and (1.3b) that for any φ ∈ C 1 ([0, T ]; R d×d

S

) and any concave or

convex g ∈ C 1 (R) d dt tr (g(φ)) = tr dφ dt g ′ (φ) = dφ dt : g ′ (φ) ∀t ∈ [0, T ]. (1.10)
Of course, a similar result holds for spatial derivatives. Furthermore, the results (1.9) and (1.10) hold true when C 1 (R) and R d×d S are replaced by C 1 (R >0 ) and R d×d S,>0 or C 1 (0, b) and R d×d S,>0,b . Finally, one can show that if g ∈ C 0,1 (R) with Lipschitz constant g Lip , then

g(φ) -g(ψ) ≤ g Lip φ -ψ ∀φ, ψ ∈ R d×d S . (1.11)
We adopt the standard notation for Sobolev spaces, e.g.

H 1 (D) := {η : D → R : D [ |η| 2 + ∇η 2 ] dx < ∞} with H 1 0 (D) being the closure of C ∞ 0 (D) for the corresponding norm • H 1 (D) . We denote the associated semi-norm as | • | H 1 (D)
. The topological dual of the Hilbert space H 1 0 (D), with pivot space L 2 (D), will be denoted by H -1 (D). Such function spaces are naturally extended when the range R is replaced by R d , R d×d and R d×d S ; e.g.

H 1 (D) becomes [H 1 (D)] d , [H 1 (D)] d×d and [H 1 (D)] d×d

S

, respectively. For ease of notation, we write the corresponding norms and semi-norms as

• H 1 (D) and | • | H 1 (D)
, respectively, as opposed to e.g. In order to analyse (P), we adopt the notation

W := [H 1 0 (D)] d , Q := L 2 (D), V := v ∈ W : D q div v dx = 0 ∀q ∈ Q , (1.12) S := [L ∞ (D)] d×d S , S >0 := [L ∞ (D)] d×d S,>0
and S >0,b := {φ ∈ S >0 : tr(φ) < b a.e. in D} .

Throughout the paper C will denote a generic positive constant independent of the regularization parameter δ and the mesh parameters h and ∆t. Finally, we recall the Poincaré inequality

(1.13) D v 2 dx ≤ C P D ∇v 2 dx ∀v ∈ W,
where C P ∈ R >0 depends only on D.

Formal free energy bound for the problem (P)

In this section we recall from Hu and Lelièvre [START_REF] Hu | New entropy estimates for the Oldroyd-B model, and related models[END_REF] the free energy structure of problem (P). Let F (u, σ) denote the free energy associated with a solution (u, p, σ) to problem (P), where we define

F (v, φ) := Re 2 D v 2 dx - ε 2Wi D b ln 1 - tr(φ) b + tr (ln(φ) + I) dx (2.1) ∀(v, φ) ∈ [L 2 (D)] d × S ⋆ with S ⋆ ⊂ S >0,b such that F (•, •) is well-defined.
Here the first term Re 2 D v 2 corresponds to the usual kinetic energy term, and the second term, which is nonnegative, is a relative entropy term. Moreover, on noting that ln is a concave function on R >0 , we observe

(2.2) F (v, φ) ≥ Re 2 D v 2 dx + ε 2Wi D tr (φ -ln(φ) -I) dx ∀(v, φ) ∈ [L 2 (D)] d × S ⋆ ,
where the right-hand side is the free energy of the Oldroyd-B model under the same no flow boundary conditions, see e.g. [START_REF] Hu | New entropy estimates for the Oldroyd-B model, and related models[END_REF] and [START_REF] Barrett | Existence and approximation of a (regularized) Oldroyd-B model[END_REF]. Clearly, diagonalization yields that the relative entropy term of this Oldroyd-B model is nonnegative. Of course, the I term in the relative entropy for FENE-P and Oldroyd-B plays no real role, and just means that the minimum relative entropy for Oldroyd-B is zero and is obtained by φ = I. Finally, we note that tr (ln(φ)) is rewritten as ln (det(φ)) in [START_REF] Hu | New entropy estimates for the Oldroyd-B model, and related models[END_REF], which once again is easily deduced from diagonalization.

Proposition 2.1. With f ∈ L 2 (0, T ; [H -1 (D)] d
) let (u, p, σ) be a sufficiently smooth solution to problem (P), (1.1a-f ), such that σ(t, •) ∈ S ⋆ for t ∈ (0, T ). Then the free energy F (u, σ) satisfies for a.a. t ∈ (0, T )

d dt F (u, σ) + (1 -ε) D ∇u 2 dx + ε 2Wi 2 D tr (A(σ)) 2 σ dx = f , u H 1 0 (D) , (2.3)
where the third term on the left-hand side is positive, via diagonalization, on recalling (1.2).

Proof. Multiplying the Navier-Stokes equation (1.1a) with u and the stress equation (1.1c) with ε 2Wi A(σ), summing and integrating over D yields, after using integrations by parts, the boundary condition (1.1f) and the incompressibility property (1.1b) in the standard way, that

D Re 2 ∂ u 2 ∂t + (1 -ε) ∇u 2 + ε Wi 1 - tr(σ) b -1 σ : ∇u dx (2.4) + ε 2Wi D ∂σ ∂t + (u • ∇)σ + A(σ) σ Wi : A(σ) dx - ε 2Wi D (∇u) σ + σ (∇u) T : A(σ) dx = f , u H 1 0 (D) .
It follows from the chain rule and (1.10) that

∂σ ∂t + (u • ∇)σ : A(σ) = ∂ ∂t + (u • ∇) -b ln 1 - tr(σ) b -tr (ln(σ)) . (2.5)
On integrating (2.5) over D, the (u • ∇) term on the right-hand side vanishes as u(t, •) ∈ V. On noting (1.2), (1.5a), (1.3b) and (1.1b), we obtain that

σ (∇u) T + (∇u) σ : A(σ) = 2 1 - tr(σ) b -1 σ : ∇u. (2.6)
Hence, on combining (2.4)-(2.6) and noting a trace property, we obtain the desired free energy equality (2.3).

For later purposes, we note the following.

Remark 2.1. The step in the above proof of testing (1.1c) with ε 2Wi A(σ) is equivalent to testing (1.1c) with -ε 2Wi σ -1 and testing the corresponding trace equation

∂ tr(σ) ∂t + (u • ∇) tr(σ) = 2∇u : σ - tr (A(σ) σ) Wi on D T (2.7) with ε 2Wi 1 -tr(σ) b -1
, and adding.

Recall that in the limit b → ∞ the FENE-P model formally converges to the Oldroyd-B model. It is thus interesting to note that when b → ∞, the free energy equality (2.3) formally converges to the corresponding free energy equality for the Oldroyd-B model, on recalling (2.2). Finally, we note the following result.

F (u(t, •), σ(t, •)) t∈(0,T ) + 1 -ε 2 DT ∇u 2 dx dt + ε 2Wi 2 DT tr (A(σ)) 2 σ dx dt (2.8) ≤ 2 F (u 0 , σ 0 ) + 1 + C P 2(1 -ε) f 2 L 2 (0,T ;H -1 (D)) .
Proof. One can bound the term f , u H 1 0 (D) in (2.3), using the Cauchy-Schwarz and Young inequalities for ν ∈ R >0 , and the Poincaré inequality (1.13), by

f , u H 1 0 (D) ≤ f H -1 (D) u H 1 (D) ≤ 1 2ν 2 f 2 H -1 (D) + ν 2 2 u 2 H 1 (D) (2.9) ≤ 1 2ν 2 f 2 H -1 (D) + ν 2 2 (1 + C P ) ∇u 2 L 2 (D)
. Combining (2.9) and (2.3) with ν 2 = (1 -ε)/(1 + C P ), and integrating in time yields the result (2.8).

3. Formal free energy bound for a regularized problem (P δ )

3.1. A regularization. Let G : s ∈ R >0 → ln s ∈ R
denote the logarithm function, whose domain of definition can be straightforwardly extended to the set of symmetric positive definite matrices using (1.6) and (1.7). We define the following concave C 1,1 (R) regularization of G based on a given parameter δ ∈ (0, 1):

G δ : s ∈ R → G(s) ∀s ≥ δ, s δ + G(δ) -1 ∀s ≤ δ ⇒ G δ (s) ≥ G(s) ∀s ∈ R >0 . (3.1)
We define also the following scalar functions

β δ (s) := (G ′ δ (s)) -1 ∀s ∈ R and β(s) := (G ′ (s)) -1 ∀s ∈ R >0 . (3.2)
Hence, we have that

β δ : s ∈ R → max{s, δ} and β : s ∈ R >0 → s. (3.3)
For later purposes, we note the following results concerning these functions. Lemma 3.1. For any φ, ψ ∈ R d×d S , η ∈ R and for any δ ∈ (0, 1) we have that

β δ (φ)G ′ δ (φ) = G ′ δ (φ)β δ (φ) = I, (3.4a) tr (ηI -G ′ δ (φ)) 2 β δ (φ) ≥ 0, (3.4b) tr (φ -G δ (φ) -I) ≥ 0, (3.4c) (φ -ψ) : [G ′ δ (ψ)] ≥ tr (G δ (φ) -G δ (ψ)) , (3.4d) -(φ -ψ) : [G ′ δ (φ) -G ′ δ (ψ)] ≥ δ 2 G ′ δ (φ) -G ′ δ (ψ) 2 . (3.4e)
In addition, if δ ∈ (0, 1 2 ] we have that

tr (φ -G δ (φ)) ≥ 1 2 φ 1 2δ [φ] - and φ : (I -G ′ δ (φ)) ≥ 1 2 φ -d. (3.5)
Proof. All the results are proved in Lemma 2.1 in [START_REF] Barrett | Existence and approximation of a (regularized) Oldroyd-B model[END_REF], except (3.4b) and this can be easily proved via diagonalization.

We introduce the following regularization of A, (1.2), for any δ ∈ (0, 1 2 ]:

A δ (φ, η) := G ′ δ 1 - η b I -G ′ δ (φ) ∀(φ, η) ∈ R d×d S × R. (3.6)
In addition to Lemma 3.1 we will also make use of the following result, which is similar to (3.5). Lemma 3.2. For any s ∈ R, b ∈ R >0 and δ ∈ (0, 1 2 ], we have that

-b G δ 1 - s b -s ≥ 1 2 [|s| -3b] + , (3.7a) G ′ δ 1 - s b -1 s ≥ [ |s| -b ] + . (3.7b)
Proof. On recalling (3.1), we first note from the concavity of G δ that

(3.8) -b G δ 1 - s b ≥ -b G δ (1) + s G ′ δ (1) = s ∀s ∈ R.
From the scalar version of (3.5), we have that

1 - s b -G δ 1 - s b ≥ 1 2 1 - s b ⇒ -b G δ 1 - s b -s ≥ b 2 1 - s b -b. (3.9) We note that b 2 1 - s b -b = 1 2 (|s| -3b) if s ≥ b, -1 2 (s + b) ≥ 1 2 (|s| -3b) if s ≤ b. (3.10) Combining (3.8)-(3.10) yields the desired result (3.7a). We now consider (3.7b). If 1 -s b ≤ δ, i.e. s ≥ b(1 -δ), then G ′ δ 1 - s b -1 s = 1 δ -1 s ≥ |s|. (3.11) If 1 -s b ≥ δ, i.e. s ≤ b(1 -δ), then G ′ δ 1 - s b -1 s = s 2 b -s ≥ 0 and s 2 b -s ≥ |s| -b. (3.12)
Combining (3.11) and (3.12) yields the desired result (3.7b).

3.2.

The regularized problem (P δ ). Using the regularizations G δ , β δ and A δ introduced above, we consider the following regularization of (P) for a given δ ∈ (0, 1 2 ]:

(P δ ) Find u δ : (t, x) ∈ [0, T ) × D → u δ (t, x) ∈ R d , p δ : (t, x) ∈ D T → p δ (t, x) ∈ R and σ δ : (t, x) ∈ [0, T ) × D → σ δ (t, x) ∈ R d×d S such that Re ∂u δ ∂t + (u δ • ∇)u δ = -∇p δ + (1 -ε)∆u δ + ε Wi div (A δ (σ δ , tr(σ δ )) β δ (σ δ )) + f on D T , (3.13a)
div u δ = 0 on D T , (3.13b)

∂σ δ ∂t + (u δ • ∇)σ δ = (∇u δ )β δ (σ δ ) + β δ (σ δ )(∇u δ ) T - A δ (σ δ , tr(σ δ )) β δ (σ δ ) Wi on D T , (3.13c
)

u δ (0, x) = u 0 (x) ∀x ∈ D, (3.13d) σ δ (0, x) = σ 0 (x) ∀x ∈ D, (3.13e) u δ = 0 on (0, T ) × ∂D. (3.13f)
3.3. Formal free energy bound for (P δ ). In this section, we extend the formal energy results (2.3) and (2.8) for (P) to problem (P δ ). We will assume throughout that

f ∈ L 2 0, T ; [H -1 (D)] d , u 0 ∈ H := {w ∈ [L 2 (D)] d : div w = 0 a.e. in D, w • n ∂D = 0 on ∂D}, (3.14) σ 0 ∈ S >0 with σ 0 min ξ 2 ≤ ξ T σ 0 (x) ξ ≤ σ 0 max ξ 2 ∀ξ ∈ R d for a.e.
x in D and tr σ 0 (x) ≤ b ⋆ for a.e. x in D;

where n ∂D is normal to ∂D, b ⋆ , σ 0 min , σ 0 max ∈ R >0 with b ⋆ < b. Let F δ (u δ , σ δ , tr(σ δ ))
denote the free energy associated with a solution (u δ , p δ , σ δ ) to problem (P δ ), where we define

F δ (v, φ, η) := Re 2 D v 2 dx - ε 2Wi D b G δ 1 - η b + tr (G δ (φ) + I) dx (3.15) ∀(v, φ, η) ∈ [L 2 (D)] d × S × L 1 (D).
Note that the second term in F δ has been regularized in comparison with F in (2.1). Similarly to (2.2), we have, on noting (3.8), the inequality (3.16)

F δ (v, φ, tr(φ)) ≥ Re 2 D v 2 dx + ε 2Wi D tr (φ -G δ (φ) -I) dx ∀(v, φ) ∈ [L 2 (D)] d × S,
where the right-hand side in (3.16) is the free energy of the corresponding regularized Oldroyd-B model, see [START_REF] Barrett | Existence and approximation of a (regularized) Oldroyd-B model[END_REF] and note (3.4c). It also follows from (3.1) and (3.14) that F δ (u 0 , σ 0 , tr(σ 0 )) ≤ F (u 0 , σ 0 ). (3.17) Proposition 3.1. Let δ ∈ (0, 1 2 ] and (u δ , p δ , σ δ ) be a sufficiently smooth solution to problem (P δ ), (3.13a-f ). Then the free energy F δ (u δ , σ δ , tr(σ δ )) satisfies for a.a. t ∈ (0, T )

d dt F δ (u δ , σ δ , tr(σ δ )) + (1 -ε) D ∇u δ 2 dx (3.18) + ε 2Wi 2 D tr (A δ (σ δ , tr(σ δ ))) 2 β δ (σ δ ) dx = f , u δ H 1 0 (D) ,
where the third term on the left-hand side is nonnegative from (3.6) and (3.4b).

Proof. Similarly to the proof of Proposition 2.1, we multiply the regularized Navier-Stokes equation (3.13a) by u δ and the regularized stress equation (3.13c) with ε 2Wi A δ (σ δ , tr(σ δ )), sum and integrate over D, use integrations by parts, the boundary condition (3.13f) and the incompressibility property (3.13b). This yields the desired result (3.18) on noting the following analogues of (2.5) and (2.6)

∂σ δ ∂t + (u δ • ∇)σ δ : A δ (σ δ , tr(σ δ )) (3.19) = ∂ ∂t + (u δ • ∇) -b G δ 1 - tr(σ δ ) b -tr(G δ (σ δ ))
and

β δ (σ δ ) (∇u δ ) T + (∇u δ ) β δ (σ δ ) : A δ (σ δ , tr(σ δ )) = 2 G ′ δ 1 - tr(σ δ ) b β δ (σ δ ) : ∇u δ . (3.20)
Here we have recalled (3.6) and (1.10) for (3.19), and (1.5a), (1.3b), (3.4a) and (3.13b) for (3.20).

Similarly to Remark 2.1, we note the following.

Remark 3.1. The step in the above proof of testing the regularized stress equation (3.13c) with

ε 2Wi A δ (σ δ , tr(σ δ )) is equivalent to testing (3.13c) with -ε 2Wi G ′ δ (σ δ
) and testing the corresponding regularized trace equation

∂ tr(σ δ ) ∂t + (u δ • ∇) tr(σ δ ) = 2∇u δ : β δ (σ δ ) - tr(A δ (σ δ , tr(σ δ )) β δ (σ δ )) Wi on D T (3.21) with ε 2Wi G ′ δ 1 -tr(σ δ ) b
, and adding.

F δ (u δ (t, •), σ δ (t, •), tr(σ δ (t, •))) + 1 -ε 2 DT ∇u δ 2 dx dt (3.22) + ε 2Wi 2 DT tr (A δ (σ δ , tr(σ δ ))) 2 β δ (σ δ ) dx dt ≤ 2 F (u 0 , σ 0 ) + 1 + C P 2(1 -ε) f 2 L 2 (0,T ;H -1 (D)) .
Proof. The proof of (3.22) follows from (3.18) in the same way as (2.8) follows from (2.3), and in addition noting (3.17).

4. Finite element approximation of (P δ ) and (P)

4.1. Finite element discretization. We now introduce a finite element discretization of the problem (P δ ), which satisfies a discrete analogue of (3.18). The time interval [0, T ) is split into intervals [t n-1 , t n ) with ∆t n = t n -t n-1 , n = 1, . . . , N T . We set ∆t := max n=1,...,NT ∆t n . We will assume throughout that the domain D is a polytope. We define a regular family of meshes {T h } h>0 with discretization parameter h > 0, which is built from partitionings of the domain D into regular open simplices so that

D = T h := NK ∪ k=1 K k with max k=1,...,NK h k ρ k ≤ C.
Here ρ k is the diameter of the largest inscribed ball contained in the simplex K k and h k is the diameter of K k , so that h = max k=1,...,NK h k . For each element K k , k = 1, . . . , N K , of the mesh T h let {P k i } d i=0 denotes its vertices, and {n k i } d i=0 the outward unit normals of the edges (d = 2) or faces (d = 3) with n k i being that of the edge/face opposite vertex P k i , i = 0, . . . , d. In addition, let {η k i (x)} d i=0 denote the barycentric coordinates of x ∈ K k with respect to the vertices {P k i } d i=0 ; that is, η k i ∈ P 1 and η k i (P k j ) = δ ij , i, j = 0, . . . , d. Here P m denote polynomials of maximal degree m in x, and δ ij the Kronecker delta notation. Finally, we introduce ∂T h := {E j } NE j=1 as the set of internal edges E j of triangles in the mesh T h when d = 2, or the set of internal faces E j of tetrahedra when d = 3.

We approximate the problem (P δ ) by the problem (P ∆t δ,h ) based on the finite element spaces W 0 h × Q 0 h × S 0 h . As is standard, we require the discrete velocity-pressure spaces W 0

h × Q 0 h ⊂ W × Q satisfy the discrete Ladyshenskaya-Babuška-Brezzi (LBB) inf-sup condition (4.1) inf q∈Q 0 h sup v∈W 0 h D q div v dx q L 2 (D) v H 1 (D) ≥ µ ⋆ > 0,
see e.g. [16, p114]. In the following, we set

W 0 h := W 2 h ⊂ W if d = 2 or W 2,- h ⊂ W if d = 2 or 3, (4.2a) Q 0 h := {q ∈ Q : q | K k ∈ P 0 k = 1, . . . , N K } ⊂ Q (4.2b) and S 0 h := {φ ∈ S : φ | K k ∈ [P 0 ] d×d S k = 1, . . . , N K } ⊂ S; (4.2c) where W 2 h := {v ∈ [C(D)] d ∩ W : v | K k ∈ [P 2 ] d k = 1, . . . , N K }, (4.3a) W 2,- h := {v ∈ [C(D)] d ∩ W : v | K k ∈ [P 1 ] d ⊕ span{ς k i } d i=0 k = 1, . . . , N K }. (4.3b)
Here, for k = 1, . . . , N K and i = 0, . . . , d

ς k i (x) = n k i d j=0,j =i η k j (x) for x ∈ K k . (4.4)
We introduce also

V 0 h := v ∈ W 0 h : D q div v dx = 0 ∀q ∈ Q 0 h , (4.5)
which approximates V. It is well-known that the choices (4.2a,b) satisfy (4.1), see e.g. [12, p221] for W 0 h = W 2 h and d = 2, and Chapter II, Sections 2.1 (d = 2) and 2.3 (d = 3) in [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations[END_REF] for W 0 h = W 2,- h . Moreover, these particular choices of S 0 h and Q 0 h have the desirable property that

(4.6) φ ∈ S 0 h ⇒ A δ (φ, tr(φ)) ∈ S 0 h and b G δ 1 - tr(φ) b + tr (G δ (φ)) ∈ Q 0 h ,
which makes it a straightforward matter to mimic the free energy inequality (3.18) at a discrete level. Since S 0 h is discontinuous, we will use the discontinuous Galerkin method to approximate the advection term (u δ • ∇)σ δ in the following. Then, for the boundary integrals, we will make use of the following definitions (see e.g. [15, p267]). Given v ∈ W 0 h , then for any φ ∈ S 0 h (or Q 0 h ) and for any point x that is in the interior of some E j ∈ ∂T h , we define the downstream and upstream values of φ at x by

φ +v (x) = lim ρ→0 + φ(x + ρ v(x)) and φ -v (x) = lim ρ→0 -φ(x + ρ v(x)); (4.7)
respectively. In addition, we denote by

[[φ]] →v (x) = φ +v (x) -φ -v (x) and {φ} v (x) = φ +v (x) + φ -v (x) 2 , (4.8)
the jump and mean value, respectively, of φ at the point x of boundary E j . From (4.7), it is clear that the values of φ +v | Ej and φ -v | Ej can change along E j ∈ ∂T h . Finally, it is easily deduced that

NE j=1 Ej |v • n|[[q 1 ]] →v q +v 2 ds = - NK k=1 ∂K k (v • n K k ) q 1 q +v 2 ds ∀v ∈ W 0 h , q 1 , q 2 ∈ Q 0 h ; (4.9)
where n ≡ n(E j ) is a unit normal to E j , whose sign is of no importance, and n K k is the outward unit normal vector of boundary ∂K k of K k . We note that similar ideas appear in upwind schemes; e.g. see Chapter IV, Section 5 in [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations[END_REF] for the Navier-Stokes equations.

A free energy preserving approximation (P ∆t

δ,h ) of (P δ ). For any source term f ∈ L 2 0, T ; [H -1 (D)] d , we define the following piecewise constant function with respect to the time variable (4.10)

f ∆t,+ (t, •) = f n (•) := 1 ∆t n t n t n-1 f (t, •) dt, t ∈ [t n-1 , t n ), n = 1, . . . , N T .
It is easily deduced that for n = 1, . . . , N T n m=1

∆t m f m r H -1 (D) ≤ t n 0 f (t, •) r H -1 (D) dt for any r ∈ [1, 2], (4.11a) and f ∆t,+ → f strongly in L 2 (0, T ; [H -1 (D)] d ) as ∆t → 0 + . (4.11b)
Throughout this section we choose u 0 h ∈ V 0 h to be the L 2 projection of u 0 onto V 0 h and σ 0 h ∈ S 0 h to be the L 2 projection of σ 0 onto S 0 h . Hence, we have that

u 0 h L 2 (D) ≤ u 0 L 2 (D) , σ 0 h | K k = 1 |K k | K k σ 0 dx, k = 1, . . . , N K , (4.12a)
where |K k | is the measure of K k ; and it immediately follows from (3.14) that

σ 0 min ξ 2 ≤ ξ T σ 0 h | K k ξ ≤ σ 0 max ξ 2 ∀ξ ∈ R d , (4.12b) tr(σ 0 h ) | K k = 1 |K k | K k tr(σ 0 ) dx ≤ tr(σ 0 ) L ∞ (K k ) ≤ b ⋆ < b. (4.12c)
We are now ready to introduce our approximation (P ∆t δ,h ) of (P δ ) for δ ∈ (0, 1 2 ]:

(P ∆t δ,h ) Setting (u 0 δ,h , σ 0 δ,h ) = (u 0 h , σ 0 h ) ∈ V 0 h × (S 0 h ∩ S >0,b ) as defined in (4.12a), then for n = 1, . . . , N T find (u n δ,h , σ n δ,h ) ∈ V 0 h × S 0 h such that for any test functions (v, φ) ∈ V 0 h × S 0 h D Re u n δ,h -u n-1 δ,h ∆t n • v + Re 2 (u n-1 δ,h • ∇)u n δ,h • v -u n δ,h • (u n-1 δ,h • ∇)v (4.13a) + (1 -ε)∇u n δ,h : ∇v + ε Wi A δ (σ n δ,h , tr(σ n δ,h )) β δ (σ n δ,h ) : ∇v dx = f n , v H 1 0 (D) , D σ n δ,h -σ n-1 δ,h ∆t n : φ -2 (∇u n δ,h ) β δ (σ n δ,h ) : φ + A δ (σ n δ,h , tr(σ n δ,h )) β δ (σ n δ,h ) : φ Wi dx (4.13b) + NE j=1 Ej u n-1 δ,h • n [[σ n δ,h ]] →u n-1 δ,h : φ +u n-1 δ,h ds = 0.
In deriving (P ∆t δ,h ), we have noted (1.5a) and that

D v • [(z • ∇)w] dx = - D w • [(z • ∇)v] dx ∀z ∈ V, ∀v, w ∈ [H 1 (D)] d , (4.14)
and we refer to [15, p267] and [START_REF] Boyaval | Free-energy-dissipative schemes for the Oldroyd-B model[END_REF] for the consistency of our approximation of the stress advection term. We note that on replacing A δ (σ n δ,h , tr(σ n δ,h )) with I -G ′ δ (σ n δ,h ) then (P ∆t δ,h ), (4.13a,b), collapses to the corresponding finite element approximation of Oldroyd-B studied in [START_REF] Barrett | Existence and approximation of a (regularized) Oldroyd-B model[END_REF], see (3.12a,b) there.

Before proving existence of a solution to (P ∆t δ,h ), we first derive a discrete analogue of the energy bound (3.18) for (P ∆t δ,h ), which uses the elementary equality

2s 1 (s 1 -s 2 ) = s 2 1 -s 2 2 + (s 1 -s 2 ) 2 ∀s 1 , s 2 ∈ R. (4.15) 4.3. Energy bound for (P ∆t δ,h ). Proposition 4.1. For n = 1, . . . , N T , a solution u n δ,h , σ n δ,h ∈ V 0 h × S 0 h to (P ∆t δ,h ), (4.13a,b), if it exists, satisfies F δ (u n δ,h , σ n δ,h , tr(σ n δ,h )) -F δ (u n-1 δ,h , σ n-1 δ,h , tr(σ n-1 δ,h )) ∆t n + Re 2∆t n D u n δ,h -u n-1 δ,h 2 dx (4.16) + (1 -ε) D ∇u n δ,h 2 dx + ε 2Wi 2 D tr A δ (σ n δ,h , tr(σ n δ,h )) 2 β δ (σ n δ,h ) dx ≤ f n , u n δ,h H 1 0 (D) ≤ (1 -ε) 2 D ∇u n δ,h 2 dx + 1 + C P 2(1 -ε) f n 2 H -1 (D) .
Proof. Similarly to (3.20), we have that

D ∇u n δ,h β δ (σ n δ,h ) : A δ (σ n δ,h , tr(σ n δ,h )) dx = D G ′ δ 1 - tr(σ n δ,h ) b β δ (σ n δ,h ) : ∇u n δ,h dx, (4.17)
where we have noted (3.6), (1.3b), (3.4a) and (4.5). Then, similarly to the proof of Proposition 3.1,

we choose v = u n δ,h ∈ V 0 h in (4.13a) and φ = ε 2Wi A δ (σ n δ,h , tr(σ n δ,h )) ∈ S 0 h in (4.13b
) and obtain, on noting (4.15), (3.4a), (4.5) and (4.17), that

f n , u n δ,h H 1 0 (D) ≥ D Re 2 u n δ,h 2 -u n-1 δ,h 2 ∆t n + u n δ,h -u n-1 δ,h 2 ∆t n + (1 -ε) ∇u n δ,h 2 dx (4.18)
+ ε 2Wi D σ n δ,h -σ n-1 δ,h ∆t n : A δ (σ n δ,h , tr(σ n δ,h )) dx + ε 2Wi 2 D tr A δ (σ n δ,h , tr(σ n δ,h )) 2 β δ (σ n δ,h ) dx + ε 2Wi NE j=1 Ej u n-1 δ,h • n [[σ n δ,h ]] →u n-1 δ,h : A δ (σ n δ,h , tr(σ n δ,h )) +u n-1 δ,h
ds.

It follows from (3.6), (3.4d) and the concavity of G δ that

σ n δ,h -σ n-1 δ,h : A δ (σ n δ,h , tr(σ n δ,h )) (4.19) ≥ tr(σ n δ,h ) -tr(σ n-1 δ,h ) G ′ δ 1 - tr(σ n δ,h ) b + tr(G δ (σ n-1 δ,h ) -tr(G δ (σ n δ,h )) ≥ b G δ 1 - tr(σ n-1 δ,h ) b + tr G δ (σ n-1 δ,h ) -b G δ 1 - tr(σ n δ,h ) b + tr G δ (σ n δ,h ) .
Similarly to (4.19), we have, on recalling (3.6), (4.7) and (4.8), that

[[σ n δ,h ]] →u n-1 δ,h : A δ (σ n δ,h , tr(σ n δ,h )) +u n-1 δ,h ≥ -[[b G δ 1 - tr(σ n δ,h ) b + tr G δ (σ n δ,h ) ]] →u n-1 δ,h . (4.20)
Finally, we note from (4.9) and as ). The second inequality in (4.16) follows immediately from (2.9) with ν 2 = (1 -ε)/(1 + C P ).

u n-1 δ,h ∈ V 0 h that for all q ∈ Q 0 h NE j=1 Ej u n-1 δ,h • n [[q]] →u n-1 δ,h ds = - NK k=1 ∂K k u n-1 δ,h • n K k q ds = - NK k=1 K k q div u n-1 δ,h dx = 0.

Existence of a solution to (P ∆t

δ,h ). Proposition 4.2. Let δ ∈ (0, 1 2 ], then, given (u n-1 δ,h , σ n-1 δ,h ) ∈ V 0 h × S 0
h and for any time step ∆t n > 0, there exists at least one solution

u n δ,h , σ n δ,h ∈ V 0 h × S 0 h to (P ∆t δ,h ), (4.13a,b).
Proof. We introduce the following inner product on the Hilbert space

V 0 h × S 0 h ((w, ψ), (v, φ)) D = D [w • v + ψ : φ] dx ∀(w, ψ), (v, φ) ∈ V 0 h × S 0 h . (4.22) Given (u n-1 δ,h , σ n-1 δ,h ) ∈ V 0 h × S 0 h , let F : V 0 h × S 0 h → V 0 h × S 0 h be such that for any (w, ψ) ∈ V 0 h × S 0 h (F (w, ψ), (v, φ)) D (4.23) := D Re w -u n-1 δ,h ∆t n • v + Re 2 (u n-1 δ,h • ∇)w • v -w • (u n-1 δ,h • ∇)v + (1 -ε)∇w : ∇v + ε Wi A δ (ψ, tr(ψ)) β δ (ψ) : ∇v + ψ -σ n-1 δ,h ∆t n : φ -2 ((∇w) β δ (ψ)) : φ + A δ (ψ, tr(ψ)) β δ (ψ) : φ Wi dx -f n , v H 1 0 (D) + NE j=1 Ej u n-1 δ,h • n [[ψ]] →u n-1 δ,h : φ +u n-1 δ,h ds ∀(v, φ) ∈ V 0 h × S 0 h .
We note that a solution (u n δ,h , σ n δ,h ) to (4.13a,b), if it exists, corresponds to a zero of F ; that is, (4.24)

F (u n δ,h , σ n δ,h ), (v, φ) D = 0 ∀(v, φ) ∈ V 0 h × S 0 h .
In addition, it is easily deduced that the mapping F is continuous. For any (w, ψ)

∈ V 0 h × S 0 h , on choosing (v, φ) = w, ε 2Wi A δ (ψ, tr(ψ))
, we obtain analogously to (4.16) that

F (w, ψ), w, ε 2Wi A δ (ψ, tr(ψ)) D (4.25) ≥ F δ (w, ψ, tr(ψ)) -F δ (u n-1 δ,h , σ n-1 δ,h , tr(σ n-1 δ,h )) ∆t n + Re 2∆t n D w -u n-1 δ,h 2 dx + 1 -ε 2 D ∇w 2 dx + ε 2Wi 2 D tr (A δ (ψ, tr(ψ))) 2 β δ (ψ) dx - 1 + C P 2(1 -ε) f n 2 H -1 (D) .
Let us now assume that for any γ ∈ R >0 , the continuous mapping F has no zero (u n δ,h , σ n δ,h ) satisfying (4.24), which lies in the ball

B γ := (v, φ) ∈ V 0 h × S 0 h : (v, φ) D ≤ γ ; (4.26) where (v, φ) D := [((v, φ), (v, φ)) D ] 1 2 = D [ v 2 + φ 2 ] dx 1 2 . (4.27)
Then for such γ, we can define the continuous mapping

G γ : B γ → B γ such that for all (v, φ) ∈ B γ G γ (v, φ) := -γ F (v, φ) F (v, φ) D . (4.28)
By the Brouwer fixed point theorem, G γ has at least one fixed point (w γ , ψ γ ) in B γ . Hence it satisfies

(4.29) (w γ , ψ γ ) D = G γ (w γ , ψ γ ) D = γ.
On noting (4.2c) and (4.29), we have that

(4.30) ψ γ 2 L ∞ (D) ≤ 1 min k∈NK |K k | D ψ γ 2 dx = µ 2 h D ψ γ 2 dx ≤ µ 2 h γ 2 ,
where

µ h := [1/(min k∈NK |K k |)]
F δ (w γ , ψ γ , tr(ψ γ )) (4.31) = Re 2 D w γ 2 dx + ε 2Wi D -b G δ 1 - tr(ψ γ ) b -tr G δ (ψ γ ) + I dx ≥ Re 2 D w γ 2 dx + ε 2Wi D tr ψ γ -G δ (ψ γ ) -I dx ≥ Re 2 D w γ 2 dx + ε 4Wi D ψ γ dx -2d|D| ≥ Re 2 D w γ 2 dx + ε 4Wi µ h γ ψ γ L ∞ (D) D ψ γ dx - εd|D| 2Wi ≥ min Re 2 , ε 4Wi µ h γ D w γ 2 + ψ γ 2 dx - εd|D| 2Wi = min Re 2 , ε 4Wi µ h γ γ 2 - εd|D| 2Wi .
Hence for all γ sufficiently large, it follows from (4.25), (4.31) and (3.4b) that (4.32)

F (w γ , ψ γ ), w γ , ε 2Wi A δ (ψ γ , tr(ψ γ )) D ≥ 0.
On the other hand as (w γ , ψ γ ) is a fixed point of G γ , we have that

F (w γ , ψ γ ), w γ , ε 2Wi A δ (ψ γ , tr(ψ γ )) D (4.33) = - F (w γ , ψ γ ) D γ D w γ 2 + ε 2Wi ψ γ : A δ (ψ γ , tr(ψ γ )) dx.
It follows from (3.6), (3.7b), (3.5) and similarly to (4.31), on noting (4.30) and (4.29), that

D w γ 2 + ε 2Wi ψ γ : A δ (ψ γ , tr(ψ γ )) dx (4.34) = D w γ 2 + ε 2Wi G ′ δ 1 - tr(ψ γ ) b -1 tr(ψ γ ) + ε 2Wi ψ γ : I -G ′ δ (ψ γ ) dx ≥ D w γ 2 dx + ε 4Wi D ψ γ dx -2d|D| ≥ min 1, ε 4Wi µ h γ γ 2 - εd|D| 2Wi .
Therefore on combining (4.33) and (4.34), we have for all γ sufficiently large that

F (w γ , ψ γ ), w γ , ε 2Wi A δ (ψ γ , tr(ψ γ )) D < 0, (4.35)
which obviously contradicts (4.32). Hence the mapping F has a zero in B γ for γ sufficiently large, and so there exists a solution (u n δ,h , σ n δ,h ) to (4.13a,b).

Theorem 4.1. For any δ ∈ (0, 1 2 ], N T ≥ 1 and any partitioning of [0, T ] into N T time steps, there exists a solution

{(u n δ,h , σ n δ,h )} NT n=1 ∈ [V 0 h × S 0 h ] NT to (P ∆t δ,h ), (4.13a,b). In addition, it follows for n = 1, . . . , N T that F δ (u n δ,h , σ n δ,h , tr(σ n δ,h )) + 1 2 n m=1 D Re u m δ,h -u m-1 δ,h 2 + (1 -ε)∆t m ∇u m δ,h 2 dx (4.36) + ε 2Wi 2 n m=1 ∆t m D tr A δ (σ m δ,h , tr(σ m δ,h )) 2 β δ (σ m δ,h ) dx ≤ F δ (u 0 h , σ 0 h , tr(σ 0 h )) + 1 + C P 2(1 -ε) n m=1 ∆t m f m 2 H -1 (D) ≤ F (u 0 h , σ 0 h ) + 1 + C P 2(1 -ε) f 2 L 2 (0,t n ;H -1 (D)) ≤ C, which yields that max n=0,...,NT D u n δ,h 2 + σ n δ,h + δ -1 [σ n δ,h ] -+ δ -1 [b -tr(σ n δ,h )] -dx ≤ C. (4.37) Moreover, for some C(h, ∆t) ∈ R >0 , but independent of δ, it follows that for k = 1, . . . , N K and n = 1, . . . , N T G ′ δ 1 - tr(σ n δ,h ) b β δ (σ n δ,h ) ≤ C(h, ∆t) on K k , (4.38a) [β δ (σ n δ,h )] -1 ≤ C(h, ∆t) 1 + G ′ δ 1 - tr(σ n δ,h ) b on K k . (4.38b)
Proof. Existence of a solution to (P ∆t δ,h ) and the first inequality in (4.36) follow immediately from Propositions 4.2 and 4.1, respectively. Similarly to (3.17), the second inequality in (4.36) is a direct consequence of (3.15), (3.1), (4.12b,c) and (4.11a). Finally, the third inequality in (4.36) follows from (4.12a-c) and (3.14).

It follows from (4.36) and (3.16) that Re 2 D u n δ,h

2 dx + ε 2Wi D tr σ n δ,h -G δ (σ n δ,h ) -I dx ≤ C, n = 1, . . . , N T . (4.39)
The first three bounds in (4.37) then follow immediately from (4.39) and (3.5). Next we note that (3.4c), (3.15) and (4.36) yield that

b D 1 - tr(σ n δ,h ) b -G δ 1 - tr(σ n δ,h ) b dx (4.40) ≤ - D b G δ 1 - tr(σ n δ,h ) b + tr G δ (σ n δ,h ) + I -b dx ≤ C.
The last bound in (4.37) is then simply obtained by using a scalar version of (3.5).

Next, we deduce from (4.36), (3.6) and (3.4b

) that for n = 1, . . . , N T 0 ≤ tr A δ (σ n δ,h , tr(σ n δ,h )) 2 β δ (σ n δ,h ) ≤ C(h, ∆t) on K k , k = 1, . . . , N K . (4.41) For any δ > 0, β δ (σ n δ,h ) ∈ R d×d S,>0
and so it follows from (1.3b), (1.4b), (3.6), (3.4a), (1.5b), (4.41), (1.8), (3.3) and (4.37) that

G ′ δ 1 - tr(σ n δ,h ) b β δ (σ n δ,h ) -I 2 (4.42) = A δ (σ n δ,h , tr(σ n δ,h )) β δ (σ n δ,h ) 2 ≤ A δ (σ n δ,h , tr(σ n δ,h )) [β δ (σ n δ,h )] 1 2 2 [β δ (σ n δ,h )] 1 2 2 = tr A δ (σ n δ,h , tr(σ n δ,h )) 2 β δ (σ n δ,h )) tr(β δ (σ n δ,h )) ≤ C(h, ∆t) β δ (σ n δ,h ) ≤ C(h, ∆t) σ n δ,h + δ ≤ C(h, ∆t) on K k , k = 1, . . . , N K , n = 1, . . . , N T .
The desired result (4.38a) follows immediately from (4.42). Similarly to (4.42), we have from (4.41), (1.8), (3.6) and (3.4a) that

C(h, ∆t) ≥ A δ (σ n δ,h , tr(σ n δ,h )) β δ (σ n δ,h ) A δ (σ n δ,h , tr(σ n δ,h )) (4.43) ≥ G ′ δ 1 - tr(σ n δ,h ) b 2 β δ (σ n δ,h ) -2G ′ δ 1 - tr(σ n δ,h ) b I + [β δ (σ n δ,h )] -1 on K k , k = 1, . . . , N K , n = 1, . . . , N T .
The desired result (4.38b) follows immediately from (4.43) and (4.38a).

4.5. Convergence of (P ∆t δ,h ) to (P ∆t h ). We now consider the corresponding direct finite element approximation of (P), i.e. (P ∆t h ) without the regularization δ:

(P ∆t h ) Given initial conditions (u 0 h , σ 0 h ) ∈ V 0 h × (S 0 h ∩ S >0,b ) as defined in (4.12a), then for n = 1, . . . , N T find (u n h , σ n h ) ∈ V 0 h × S 0 h such that for any test functions (v, φ) ∈ V 0 h × S 0 h D Re u n h -u n-1 h ∆t n • v + Re 2 (u n-1 h • ∇)u n h • v -u n h • (u n-1 h • ∇)v (4.44a) + (1 -ε)∇u n h : ∇v + ε Wi A(σ n h ) σ n h : ∇v dx = f n , v H 1 0 (D) , D σ n h -σ n-1 h ∆t n : φ -2 ((∇u n h ) σ n h ) : φ + A(σ n h ) σ n h : φ Wi dx (4.44b) + NE j=1 Ej u n-1 h • n [[σ n h ]] →u n-1 h : φ +u n-1 h ds = 0.
We note that (4.44a,b) and

F (u n h , σ n h ) are only well-defined if σ n h ∈ S 0 h ∩ S >0,b
. We also note that on replacing A(σ n h ) with I -(σ n h ) -1 then (P ∆t h ), (4.44a,b), collapses to the corresponding finite element approximation of Oldroyd-B studied in [START_REF] Barrett | Existence and approximation of a (regularized) Oldroyd-B model[END_REF], see (

{(u n h , σ n h )} NT n=1 ∈ [V 0 h × S 0 h ]
NT such that for the subsequence

u n δ,h → u n h , σ n δ,h → σ n h as δ → 0 + , for n = 1, . . . , N T . (4.45) In addition, for all t ∈ [0, T ] n = 1, . . . , N T , σ n h | K k ∈ R d×d S,>0,b , k = 1, . . . , N K ,. Moreover, {(u n h , σ n h )} NT n=1 ∈ [V 0 h × S 0 h ] NT solves (P ∆t h ), (4.44a,b), and for n = 1, . . . , N T F (u n h , σ n h ) -F (u n-1 h , σ n-1 h ) ∆t n + Re 2∆t n D u n h -u n-1 h 2 dx + (1 -ε) D ∇u n h 2 dx (4.46) + ε 2Wi 2 D tr (A(σ n h )) 2 σ n h dx ≤ 1 2 (1 -ε) D ∇u n h 2 dx + 1 + C P 2(1 -ε) f n 2 H -1 (D) .
Proof

β δ (σ n h ) → σ n h , β δ (σ n δ,h ) → σ n h as δ → 0 + . (4.47) If tr(σ n h )| K k = b on some simplex K k , then for the subsequence of (4.45) we have that tr(σ n δ,h )| K k → b as δ → 0 + . (4.48)
In addition, it follows from (1.8), (3.3) and (4.48) for δ sufficiently small that

β δ (σ n δ,h ) ≥ 1 √ d tr β δ (σ n δ,h ) ≥ 1 √ d tr σ n δ,h ≥ b 2 √ d on K k . (4.49)
Hence, (4.38a) and (4.49) yield for the subsequence of (4.45) for all δ sufficiently small that

G ′ δ 1 - tr(σ n δ,h ) b ≤ C(h, ∆t) on K k , ( 4
G δ 1 - tr(σ n δ,h ) b → ln 1 - tr(σ n h ) b , G ′ δ 1 - tr(σ n δ,h ) b → 1 - tr(σ n h ) b -1 (4.51) as δ → 0 + .
The results (4.38b) and (4.51) yield for δ sufficiently small that

[β δ (σ n δ,h )] -1 ≤ C(h, ∆t) on K k , k = 1, . . . , N K . (4.52)
Furthermore, it follows from (4.52), (4.47), (3.2), (3.3) and as

[β δ (σ n δ,h )] -1 β δ (σ n δ,h ) = I that the following subsequence result [β δ (σ n δ,h )] -1 = G ′ δ (σ n δ,h ) → [σ n h ] -1 as δ → 0 + (4.53)
holds, and so

σ n h | K k ∈ R d×d S,>0,b , k = 1, . . . ,
∈ S 0 h Ej u n-1 δ,h • n [[σ n δ,h ]] →u n-1 δ,h : φ +u n-1 δ,h ds → Ej u n-1 h • n [[σ n h ]] →u n-1 h : φ +u n-1 h ds as δ → 0 + . (4.55)
Hence using (4.45), (4.47), (4.51) and (4.55), we can pass to the limit δ → 0 + for the subsequence in

(P ∆t δ,h ), (4.13a,b), to show that {(u n h , σ n h )} NT n=1 ∈ [V 0 h × S 0 h ] NT solves (P ∆t h ), (4.44a
,b). Similarly, using (4.45), (4.47), (4.51), (4.53) and (4.54), and noting (3.15) and (2.1), we can pass to the limit δ → 0 + in (4.16) to obtain (4.46).

Fene-p model with stress diffusion

5.1. Model (P α ), (P) with stress diffusion. In this section we consider the following modified version of (P), (1.1a-f), with a stress diffusion term for a given constant α ∈ R >0 :

(P α ) Find u α : (t, x) ∈ [0, T ) × D → u α (t, x) ∈ R d , p α : (t, x) ∈ D T → p α (t, x) ∈ R and σ α : (t, x) ∈ [0, T ) × D → σ α (t, x) ∈ R d×d S,>0,b such that Re ∂u α ∂t + (u α • ∇)u α = -∇p α + (1 -ε)∆u α + ε Wi div (A(σ α ) σ α ) + f on D T , (5.1a) div u α = 0 on D T , (5.1b) ∂σ α ∂t + (u α • ∇)σ α = (∇u α )σ α + σ α (∇u α ) T - A(σ α ) σ α Wi + α∆σ α on D T , (5.1c) u α (0, x) = u 0 (x) ∀x ∈ D, (5.1d) σ α (0, x) = σ 0 (x) ∀x ∈ D, (5.1e) u α = 0 on (0, T ) × ∂D, (5.1f) (n ∂D • ∇)σ α = 0 on (0, T ) × ∂D. (5.1g)
Hence problem (P α ) is the same as (P), but with the added diffusion term α∆σ α for the stress equation (5.1c), and the associated Neumann boundary condition (5.1g). Similarly to (P δ ), (3.13a-f), we introduce a regularization of (P α,δ ) of (P α ) mimicking the free energy structure of (P α ). Moreover, we need to be able to construct a finite element approximation of (P α,δ ) that satisfies a discrete analogue of this free energy structure. Apart from the obvious addition of the stress diffusion term, there are three other distinct differences. First, one has to deal with the advective term in the stress equation differently, see (5.11) below, to the approach used in (3.19) for the stability bound, as the approach there cannot be mimicked at a discrete level using continuous piecewise linear functions to approximate σ α,δ , the regularization of σ α . Note that one cannot use S 0 h with the desirable property (4.6) to approximate σ α,δ due to the additional stress diffusion term. Second, as a consequence of this stress advective term, one has to introduce another regularization of tr(σ α ), ̺ α,δ , as well as the obvious candidate tr(β δ (σ α,δ )), and solve for this directly, see Remark 5.1 below. Third, it is desirable for the convergence analysis, as δ → 0, to have a uniform L 2 (D T ) bound on the extra stress term A δ (σ α,δ , ̺ α,δ ) β δ (σ α,δ ) in the Navier-Stokes equation (5.5a) below. To achieve this we replace A δ (σ α,δ , ̺ α,δ ) β δ (σ α,δ ) there by κ δ (σ α,δ , ̺ α,δ ) A δ (σ α,δ , ̺ α,δ ) β δ (σ α,δ ), where we define, for δ ∈ (0, min{ 1 2 , b}],

κ δ (φ, η) := β b δ (η) tr(β δ (φ)) 1 2 ∀(φ, η) ∈ R d×d S × R (5.2) with β b δ : s ∈ R → min{β δ (s), b}. (5.
3) It follows from (5.2), (5.3), (1.3b) and (1.4b) that

κ δ (φ, η) A δ (φ, η) β δ (φ) 2 ≤ κ δ (φ, η) [β δ (φ)] 1 2 2 A δ (φ, η) [β δ (φ)] 1 2 2 (5.4) ≤ b tr (A δ (φ, η)) 2 β δ (φ) ∀(φ, η) ∈ R d×d S × R. Hence a uniform L 2 (D T ) bound on κ δ (σ α,δ , ̺ α,δ ) A δ (σ α,δ , ̺ α,δ ) β δ (σ α,δ ) follows from a uniform L 1 (D T ) on tr (A δ (σ α,δ , ̺ α,δ )) 2 β δ (σ α,δ
) , which will follow from the free energy bound. Although

̺ α,δ = tr(β δ (σ α,δ
)), we will show, in the limit δ → 0, that

β δ (σ α,δ ) → σ α and ̺ α,δ → tr(σ α ) with σ α (•, •) ∈ R d×d S,>0,b
, and hence implying that κ δ (σ α,δ , ̺ α,δ ) → 1. In order to maintain the free energy bound we need to include κ δ (σ α,δ , ̺ α,δ ) on the right-hand sides of (5.5c,d) below.

Therefore, we consider the following regularization of (P α ) for a given δ ∈ (0, min{ 1 2 , b}]:

(P α,δ ) Find u α,δ : (t, x) ∈ [0, T ) × D → u α,δ (t, x) ∈ R d , p α,δ : (t, x) ∈ D T → p α,δ (t, x) ∈ R, σ α,δ : (t, x) ∈ [0, T ) × D → σ α,δ (t, x) ∈ R d×d S and ̺ α,δ : (t, x) ∈ [0, T ) × D → ̺ α,δ (t, x) ∈ R such that Re ∂u α,δ ∂t + (u α,δ • ∇)u α,δ = -∇p α,δ + (1 -ε)∆u α,δ (5.5a) + ε Wi div (κ δ (σ α,δ , ̺ α,δ )A δ (σ α,δ , ̺ α,δ ) β δ (σ α,δ )) + f on D T , div u α,δ = 0 on D T , (5.5b) ∂σ α,δ ∂t + (u α,δ • ∇)β δ (σ α,δ ) = κ δ (σ α,δ , ̺ α,δ ) (∇u α,δ )β δ (σ α,δ ) + β δ (σ α,δ )(∇u α,δ ) T (5.5c) - A δ (σ α,δ , ̺ α,δ ) β δ (σ α,δ ) Wi + α∆σ α,δ on D T , ∂̺ α,δ ∂t -b (u α,δ • ∇)β δ (1 - ̺ α,δ b ) = 2κ δ (σ α,δ , ̺ α,δ ) ∇u α,δ : β δ (σ α,δ ) (5.5d) - tr(A δ (σ α,δ , ̺ α,δ ) β δ (σ α,δ )) Wi + α∆̺ α,δ on D T , u α,δ (0, x) = u 0 (x) ∀x ∈ D, (5.5e) σ α,δ (0, x) = σ 0 (x), ̺ α,δ (0, x) = tr(σ 0 (x)) ∀x ∈ D, (5.5f) u α,δ = 0 on (0, T ) × ∂D, (5.5g) (n ∂D • ∇)σ α,δ = 0, (n ∂D • ∇)̺ α,δ = 0 on (0, T ) × ∂D. (5.5h) We note from (3.3) that -b ∇β δ 1 - ̺ α,δ b = ∇̺ α,δ if 1 - ̺ α,δ b ≥ δ, i.e. ̺ α,δ ≤ b(1 -δ).
We remark again, due the required regularization of the advective terms in (5.5c,d), that ̺ α,δ = tr(σ α,δ ). However, we note that on taking the trace of (5.5c), subtracting (5.5d) and integrating over D, yields, on noting (5.5b,f-h) and (1.3b), that for all t ∈ [0, T ]

D [tr(σ α,δ (t, •)) -̺ α,δ (t, •)] dx = D [tr(σ α,δ (0, •)) -̺ α,δ (0, •)] dx = 0.
(5.6) 5.2. Formal free energy bound for (P α,δ ). First, similarly to (3.1), we introduce, for δ ∈ (0, 1), the concave C 1,1 (R >0 ) function

H δ : s ∈ R >0 → G(s) ∀s ∈ (0, δ -1 ], δ s + G(δ -1 ) -1 ∀s ≥ δ -1 ⇒ H ′ δ (G ′ δ (s)) = β δ (s) ∀s ∈ R. (5.7)
We have the following analogue of Proposition 3.1. 1 2 , b}] and (u α,δ , p α,δ , σ α,δ , ̺ α,δ ) be a sufficiently smooth solution to problem (P α,δ ), (5.5a-h). Then the free energy

Proposition 5.1. Let α ∈ R >0 , δ ∈ (0, min{
F δ (u α,δ , σ α,δ , ̺ α,δ ) satisfies for a.a. t ∈ (0, T ) d dt F δ (u α,δ , σ α,δ , ̺ α,δ ) + (1 -ε) D ∇u α,δ 2 dx (5.8) + αεδ 2 2Wi D ∇G ′ δ (σ α,δ ) 2 + b ∇G ′ δ 1 - ̺ α,δ b 2 dx + ε 2Wi 2 D tr (A δ (σ α,δ , ̺ α,δ )) 2 β δ (σ α,δ ) dx ≤ f , u α,δ H 1 0 (D) .
Proof. Similarly to the proof of Proposition 3.18, on noting Remark 3.1, we multiply the Navier-Stokes equation (5.5a) with u α,δ and the stress equation (5.5c) with

-ε 2Wi G ′ δ (σ α,δ ), the trace equation (5.5d) with ε 2Wi G ′ δ 1 - ̺ α,δ b
, sum and integrate over D. This yields, after performing integration by parts and noting (5.5b,g,h), (1.5a), (3.4a) and (3.6), that

D Re 2 ∂ u α,δ 2 ∂t + (1 -ε) ∇u α,δ 2 dx (5.9) - ε 2Wi D ∂σ α,δ ∂t + (u α,δ • ∇)β δ (σ α,δ ) : G ′ δ (σ α,δ ) dx + ε 2Wi D ∂̺ α,δ ∂t -b (u α,δ • ∇)β δ 1 - ̺ α,δ b G ′ δ 1 - ̺ α,δ b dx - αε 2Wi D ∇σ α,δ :: ∇G ′ δ (σ α,δ ) -∇̺ α,δ • ∇G ′ δ 1 - ̺ α,δ b dx + ε 2Wi 2 D tr (A δ (σ α,δ , ̺ α,δ )) 2 β δ (σ α,δ ) dx = f , u α,δ H 1 0 (D) . Using (1.10), we have that ∂σ α,δ ∂t : G ′ δ (σ α,δ ) = ∂ ∂t tr (G δ (σ α,δ )) and ∂̺ α,δ ∂t G ′ δ 1 - ̺ α,δ b = -b ∂ ∂t G δ 1 - ̺ α,δ b .
(5.10)

As β δ (σ α,δ ) ≡ H ′ δ (G ′ δ (σ α,δ
)), on recalling (5.7), we have, on noting u α,δ ∈ V and the spatial version of (1.10), that

- D (u α,δ • ∇)β δ (σ α,δ ) : G ′ δ (σ α,δ ) dx = D β δ (σ α,δ ) : (u α,δ • ∇)G ′ δ (σ α,δ ) dx (5.11) = D (u α,δ • ∇) tr (H δ (G ′ δ (σ α,δ ))) dx = 0.
Similarly to (5.11), we have that

- D (u α,δ • ∇)β δ 1 - ̺ α,δ b G ′ δ 1 - ̺ α,δ b dx = D (u α,δ • ∇)H δ G ′ δ 1 - ̺ α,δ b dx = 0.
(5.12)

Similarly to (3.4e), we have that -∇σ α,δ :: ) I and integrate over D, then we obtain, on noting u α,δ ∈ V and (5.12), that

∇G ′ δ (σ α,δ ) ≥ δ 2 ∇G ′ δ (σ α,δ ) 2 a.e. in D T , (5.13a) ∇̺ α,δ • ∇G ′ δ 1 - ̺ α,δ b = -b∇ 1 - ̺ α,δ b • ∇G ′ δ 1 - ̺ α,δ b (5.13b) ≥ b δ 2 ∇G ′ δ 1 - ̺ α,
- D (u α,δ • ∇)β δ (σ α,δ ) : G ′ δ 1 - tr(σ α,δ ) b Idx = D tr(β δ (σ α,δ )) (u α,δ • ∇)G ′ δ 1 - tr(σ α,δ ) b dx = -b D 1 - tr(β δ (σ α,δ )) b (u α,δ • ∇)G ′ δ 1 - tr(σ α,δ ) b dx = -b D β δ 1 - tr(σ α,δ ) b (u α,δ • ∇)G ′ δ 1 - tr(σ α,δ ) b dx = 0.
Hence, the need for the new variable ̺ α,δ in order to mimic the free energy structure of (P α ).

The following Corollary follows from (5.8) on noting the proof of Corollary 3.1.

Corollary 5.1. Under the assumptions of Proposition 5.1 it follows that sup t∈(0,T )

F δ (u α,δ (t, •), σ α,δ (t, •), ̺ α,δ (t, •)) + 1 -ε 2 DT ∇u α,δ 2 dx dt (5.14) + αεδ 2 2Wi DT ∇G ′ δ (σ α,δ ) 2 + b ∇G ′ δ 1 - ̺ α,δ b 2 dx dt + ε 2Wi 2 DT tr (A δ (σ α,δ , ̺ α,δ )) 2 β δ (σ α,δ ) dx dt ≤ 2 F (u 0 , σ 0 ) + 1 + C P 2(1 -ε) f 2 L 2 (0,T ;H -1 (D)) .
6. Finite element approximation of (P α,δ )

6.1. Finite element discretization. We now introduce a conforming finite element discretization of (P α,δ ), (5.5a-h), which satisfies a discrete analogue of (5.8). As noted in Section 5, we cannot use S 0 h with the desirable property (4.6) to approximate σ α,δ , as we now have the added diffusion term. In the following, we choose

W 1 h := W 2 h ⊂ W or W 1,+ h ⊂ W, (6.1a) Q 1 h = {q ∈ C(D) : q | K k ∈ P 1 k = 1, . . . , N K } ⊂ Q, (6.1b) S 1 h = {φ ∈ [C(D)] d×d S : φ | K k ∈ [P 1 ] d×d S k = 1, . . . , N K } ⊂ S (6.1c) and V 1 h = v ∈ W 1 h : D q div v dx = 0 ∀q ∈ Q 1 h ; (6.1d)
where W 2 h is defined as in (4.3a) and, on recalling the barycentric coordinate notation used in (4.4),

W 1,+ h :=    v ∈ [C(D)] d ∩ W : v | K k ∈ P 1 ⊕ span d i=0 η k i d k = 1, . . . , N K    . (6.2)
The velocity-pressure choice, W 2 h × Q 1 h , is the lowest order Taylor-Hood element. It satisfies (4.1) with W 0 h and Q 0 h replaced by W 2 h and Q 1 h , respectively, provided, in addition to {T h } h>0 being a regular family of meshes, that each simplex has at least one vertex in D, see p177 in Girault and Raviart [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations[END_REF] in the case d = 2 and Boffi [START_REF] Boffi | Three-dimensional finite element methods for the Stokes problem[END_REF] in the case d = 3. Of course, this is a very mild restriction on {T h } h>0 . The velocity-pressure choice, W [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations[END_REF] in the case d = 2 and Section 4.2.4 in Ern and Guermond [START_REF] Ern | Theory and Practice of Finite Elements[END_REF] in the case d = 3. Hence for both choices of W 1 h , it follows that for all v ∈ V there exists a sequence

{v h } h>0 , with v h ∈ V 1 h , such that lim h→0+ v -v h H 1 (D) = 0 . (6.3)
We recall the well-known local inverse inequality for Q

1 h q L ∞ (K k ) ≤ C |K k | -1 K k |q| dx ∀q ∈ Q 1 h , k = 1, . . . , N K (6.4) ⇒ χ L ∞ (K k ) ≤ C |K k | -1 K k χ dx ∀χ ∈ S 1 h , k = 1, . . . , N K .
We recall a similar well-known local inverse inequality for S 1

h and V 1 h ∇φ L 2 (K k ) ≤ C h -1 k φ L 2 (K k ) ∀φ ∈ S 1 h , k = 1, . . . , N K , (6.5a) ∇v L 2 (K k ) ≤ C h -1 k v L 2 (K k ) ∀v ∈ V 1 h , k = 1, . . . , N K . (6.5b)
We introduce the interpolation operator π h : C(D) → Q 1 h , and extended naturally to

π h : [C(D)] d×d S → S 1 h , such that for all η ∈ C(D) and φ ∈ [C(D)] d×d S π h η(P p ) = η(P p ) and π h φ(P p ) = φ(P p ) p = 1, . . . , N P , (6.6)
where {P p } NP p=1 are the vertices of

T h . As φ ∈ S 1 h and q ∈ Q 1 h do not imply that G ′ δ (φ) ∈ S 1 h and G ′ δ 1 -q b ∈ Q 1 h
, we have to test the finite element approximation of (5.5c,d) with

-ε 2Wi π h [G ′ δ (σ n α,δ,h )] ∈ S 1 h and ε 2Wi π h G ′ δ 1 - ̺ n α,δ,h b ∈ Q 1 h
, respectively, where σ n α,δ,h ∈ S 1 h and ̺ n α,δ,h ∈ Q 1 h are our finite element approximations to σ α,δ and ̺ α,δ at time level t n . Therefore the finite element approximation of (5.5c,d) have to be constructed to mimic the results (5.9)-(5.13a,b), when tested with

-ε 2Wi π h [G ′ δ (σ n α,δ,h )] ∈ S 1 h and ε 2Wi π h G ′ δ 1 - ̺ n α,δ,h b ∈ Q 1 h , respectively.
In order to mimic (5.10) and the (P α,δ ) analogue of (3.20), we need to use numerical integration (vertex sampling). We note the following results. As the basis functions associated with Q 1 h and S 1 h are nonnegative and sum to unity everywhere, we have, on noting (1.5b), for k = 1, . . . , N K that

π h [φ ψ] ≤ π h [ φ ψ ] (6.7a) ≤ [π h [ φ r1 ] ] 1 r 1 [π h [ ψ r2 ] ] 1 r 2 on K k , ∀φ, ψ ∈ [C(K k )] d×d S , π h φ 2 ≤ π h [ φ 2 ] on K k , ∀φ ∈ [C(K k )] d×d S , (6.7b) where r 1 , r 2 ∈ (1, ∞) satisfy r -1 1 + r -1 2 = 1.
In addition, we have for k = 1, . . . , N K that

K k χ 2 dx ≤ K k π h [ χ 2 ] dx ≤ C K k χ 2 dx ∀χ ∈ S 1 h . (6.8)
The first inequality in (6.8) follows immediately from (6.7b), and the second from applying (6.4) and a Cauchy-Schwarz inequality. Of course, scalar versions of (6.7a,b) and (6.8) hold with [C(K k )] d×d S and S 1 h replaced by C(K k ) and Q 1 h , respectively. Furthermore, for later use, we recall the following well-known results concerning the interpolant π h for k = 1, . . . , N K :

(I -π h )q W 1,∞ (K k ) ≤ C h |q| W 2,∞ (K k ) ∀q ∈ W 2,∞ (K k ), (6.9a) (I -π h )[q 1 q 2 ] L 1 (K k ) ≤ C h 2 k ∇q 1 L 2 (K k ) ∇q 2 L 2 (K k ) (6.9b) ≤ C h k q 1 L 2 (K k ) ∇q 2 L 2 (K k ) ∀q 1 , q 2 ∈ Q 1 h .
In order to mimic (5.11) and (5.12), we have to carefully construct our finite element approximation of the advective terms in (5.5c,d). Our construction is a non-trivial extension of an approach that has been used in the finite element approximation of fourth-order degenerate nonlinear parabolic equations, such as the thin film equation; see e.g. Grün and Rumpf [START_REF] Grün | Nonnegativity preserving numerical schemes for the thin film equation[END_REF] and Barrett and Nürnberg [START_REF] Barrett | Convergence of a finite-element approximation of surfactant spreading on a thin film in the presence of van der Waals forces[END_REF]. Let {e i } d i=1 be the orthonormal vectors in R d , such that the j th component of e i is δ ij , i, j = 1, . . . , d. Let K be the standard open reference simplex in R d with vertices { P i } d i=0 , where P 0 is the origin and P i = e i , i = 1, . . . , d. Given a simplex K k ∈ T h with vertices {P k i } d i=0 , then there exists a non-singular matrix B k such that the linear mapping (6.10)

B k : x ∈ R d → P k 0 + B k x ∈ R d maps vertex P i to vertex P k i , i = 0, . . . , d. Hence B k maps K to K k . For all η ∈ Q 1 h and K k ∈ T h , we define η( x) := η(B k ( x)) ∀ x ∈ K ⇒ ∇η(B k ( x)) = (B T k ) -1 ∇ η( x) ∀ x ∈ K , (6.11)
where for all x ∈ K

[ ∇ η( x)] j = ∂ ∂ x j η( x) = η( P j ) -η( P 0 ) = η(P k j ) -η(P k 0 ) j = 1, . . . , d. (6.12)
Such notation is easily extended to φ ∈ S 1 h . Given q ∈ Q 1 h and K k ∈ T h , then first, for j = 1, . . . , d, we find the unique Λ k δ,j ( q) ∈ R, which is continuous on q, such that

Λ k δ,j ( q) ∂ ∂ x j π h [G ′ δ ( q)] = ∂ ∂ x j π h [H δ (G ′ δ ( q))] on K, (6.13)
where (

π h η)( x) ≡ (π h η)(B k x) for all x ∈ K and η ∈ C(K k ). We set Λ k δ,j ( q) :=      H δ (G ′ δ (q(P k j ))) -H δ (G ′ δ (q(P k 0 ))) G ′ δ (q(P k j )) -G ′ δ (q(P k 0 )) if β δ (q(P k j )) = β δ (q(P k 0 )), β δ (q(P k j )) = β δ (q(P k 0 )) if β δ (q(P k j )) = β δ (q(P k 0 )), (6.14)
where we have noted that G ′ δ (β δ (s)) = G ′ δ (s) for all s ∈ R and G ′ δ (•) is strictly decreasing on [δ, ∞). Clearly, Λ k δ,j ( q) ∈ R, j = 1, . . . , d, satisfies (6.13) and depends continuously on q | K k . Next, we extend the construction (6.14) for a given φ ∈ S 1 h and K k ∈ T h , to find for j = 1, . . . , d the unique Λ k δ,j ( φ) ∈ R d×d

S

, which is continuous on φ, such that

Λ k δ,j ( φ) : ∂ ∂ x j π h [G ′ δ ( φ)] = ∂ ∂ x j π h [tr(H δ (G ′ δ ( φ)))] on K. (6.15)
To construct Λ k δ,j ( φ) satisfying (6.15), we note the following. We have from (5.7) and (1.9) that

β δ (φ(P k j )) : (G ′ δ (φ(P k j )) -G ′ δ (φ(P k 0 ))) ≤ tr(H δ (G ′ δ (φ(P k j )) -H δ (G ′ δ (φ(P k 0 ))) (6.16) ≤ β δ (φ(P k 0 )) : (G ′ δ (φ(P k j )) -G ′ δ (φ(P k 0 ))). Since G ′ δ (β δ (s)) = G ′ δ (s) for all s ∈ R, it follows from (3.4e) that -(β δ (φ(P k j )) -β δ (φ(P k 0 ))) : (G ′ δ (φ(P k j )) -G ′ δ (φ(P k 0 ))) (6.17) = -(β δ (φ(P k j )) -β δ (φ(P k 0 ))) : (G ′ δ (β δ (φ(P k j ))) -G ′ δ (β δ (φ(P k 0 )))) ≥ δ 2 G ′ δ (β δ (φ(P k j ))) -G ′ δ (β δ (φ(P k 0 ))) 2 .
Therefore the left-hand side of (6.17) is zero if and only if ), the range of β δ (•). Hence, on noting (6.12), (6.16), (6.17) and (1.3b), we have that

G ′ δ (β δ (φ(P k j ))) = G ′ δ (β δ (φ(P k 0 ))), which is equivalent to β δ (φ(P k j )) = β δ (φ(P k 0 )) as G ′ δ (•) is invertible on [δ, ∞
Λ k δ,j ( φ) := (1 -λ k δ,j )β δ (φ(P k j )) + λ k δ,j β δ (φ(P k 0 )) (6.18a) if (β δ (φ(P k j )) -β δ (φ(P k 0 ))) : (G ′ δ (φ(P k j )) -G ′ δ (φ(P k 0 ))) = 0 , Λ k δ,j ( φ) := β δ (φ(P k j )) = β δ (φ(P k 0 )) (6.18b) if (β δ (φ(P k j )) -β δ (φ(P k 0 ))) : (G ′ δ (φ(P k j )) -G ′ δ (φ(P k 0 )
)) = 0 satisfies (6.15) for j = 1, . . . , d; where λ k δ,j ∈ [0, 1] is defined as

λ k δ,j := tr(H δ (G ′ δ (φ(P k j )) -H δ (G ′ δ (φ(P k 0 ))) -β δ (φ(P k j )) : (G ′ δ (φ(P k j )) -G ′ δ (φ(P k 0 ))) (β δ (φ(P k 0 )) -β δ (φ(P k j ))) : (G ′ δ (φ(P k j )) -G ′ δ (φ(P k 0 ))) . Furthermore, Λ k δ,j ( φ) ∈ R d×d S , j = 1, . . . , d, depends continuously on φ | K k . Therefore given q ∈ Q 1 h and φ ∈ S 1 h , we introduce, for m, p = 1, . . . , d, Λ δ,m,p (q) = d j=1 [(B T k ) -1 ] mj Λ k δ,j ( q) [B T k ] jp ∈ R, on K k , k = 1, . . . , N K , (6.19a) Λ δ,m,p (φ) = d j=1 [(B T k ) -1 ] mj Λ k δ,j ( φ) [B T k ] jp ∈ R d×d S on K k , k = 1, . . . , N K . (6.19b)
It follows from (6.19a,b), (6.13), (6.15) and (6.11) that Λ δ,m,p (q) ≈ β δ (q) δ mp , Λ δ,m,p (φ) ≈ β δ (φ) δ mp m, p = 1, . . . , d ; (6.20) and for m = 1, . . . , d

d p=1 Λ δ,m,p (q) ∂ ∂x p π h [G ′ δ (q)] = ∂ ∂x m π h [H δ (G ′ δ (q))] on K k , k = 1, . . . , N K , (6.21a) d p=1 Λ δ,m,p (φ) : ∂ ∂x p π h [G ′ δ (φ)] = ∂ ∂x m π h [tr(H δ (G ′ δ (φ)))] on K k , k = 1, . . . , N K . (6.21b)
For a more precise version of (6.20), see Lemma 6.4 below. Of course, for (6.19a) and (6.21a) we can adopt the more compact notation on

K k , k = 1, . . . , N K , Λ δ (q) = (B T k ) -1 Λ k δ ( q) B T k ≈ β δ (q) I ⇒ Λ δ (q)∇π h [G ′ δ (q)] = ∇π h [H δ (G ′ δ (q))], (6.22) where Λ k δ ( q) ∈ R d×d S is diagonal with [Λ k δ ( q)] jj = Λ k δ,j (
q), j = 1, . . . , d, so that Λ δ (q) ∈ R d×d S with [Λ δ (q)] mp = Λ δ,m,p (q), m, p = 1, . . . , d.

Finally, as the partitioning T h consists of regular simplices, we have that

(B T k ) -1 B T k ≤ C, k = 1, . . . , N K . (6.23)
Hence, it follows from (6.19a,b), (6.23) and (6.18a,b) that for k

= 1, . . . , N K Λ δ,m,p (q) L ∞ (K k ) ≤ C π h [β δ (q)] L ∞ (K k ) ∀q ∈ Q 1 h , (6.24a) Λ δ,m,p (φ) L ∞ (K k ) ≤ C π h [β δ (φ)] L ∞ (K k ) ∀φ ∈ S 1 h . (6.24b)
In order to mimic (5.13a,b), we shall assume from now on that the family of meshes, {T h } h>0 , for the polytope D consists of non-obtuse simplices only, i.e. all dihedral angles of any simplex in T h are less than or equal to π 2 . Let K k have vertices {P k j } d j=0 , and let η k j (x) be the basis functions on

K k associated with Q 1 h and S 1 h , i.e. η k j | K k ∈ P 1 and η k j (P k i ) = δ ij , i, j = 0, . . . , d. As K k is non-obtuse it follows that ∇η k i • ∇η k j ≤ 0 on K k , i, j = 0, . . . , d, with i = j. (6.25)
We then have the following result. Lemma 6.1. Let g ∈ C 0,1 (R) be monotonically increasing with Lipschitz constant g Lip . As T h consists of only non-obtuse simplices, then we have for all q ∈ Q 1

h , φ ∈ S 1 h that g Lip ∇π h [g(q)] • ∇q ≥ ∇π h [g(q)] 2 and g Lip ∇π h [g(φ)] :: ∇φ ≥ ∇π h [g(φ)] 2 (6.26) on K k , k = 1, . . . , N K .
Proof. See the proof of Lemma 5.1 in [START_REF] Barrett | Existence and approximation of a (regularized) Oldroyd-B model[END_REF].

Of course, the construction of a non-obtuse mesh in the case d = 3 is not straightforward for a general polytope D. However, we stress that our numerical method (P ∆t α,δ,h ), see (6.34a-c) below, does not require this constraint. It is only required to show that (P ∆t α,δ,h ) mimics the free energy structure of (P α,δ ).

6.2.

A free energy preserving approximation, (P ∆t α,δ,h ), of (P α,δ ). In addition to the assumptions on the finite element discretization stated in subsection 6.1, and our definition of ∆t in subsection 4.1, we shall assume for the convergence analysis, see Section 7, that there exists a

C ∈ R >0 such that ∆t n ≤ C ∆t n-1 , n = 2, . . . , N, as ∆t → 0 + . (6.27)
We note that this constraint is not required for the results in this section, in particular Theorem 6.1.

With ∆t 1 and C as above, let ∆t 0 ∈ R >0 be such that ∆t 1 ≤ C∆t 0 . Given initial data satisfying (3.14), we choose u 0 h ∈ V 1 h and σ 0 h ∈ S 1 h throughout the rest of this paper such that

D u 0 h • v + ∆t 0 ∇u 0 h : ∇v dx = D u 0 • v dx ∀v ∈ V 1 h , (6.28a) D π h [σ 0 h : χ] + ∆t 0 ∇σ 0 h :: ∇χ dx = D σ 0 : χ dx ∀χ ∈ S 1 h . (6.28b)
It follows from (6.28a,b), (6.8) and (3.14) that

D u 0 h 2 + σ 0 h 2 + ∆t 0 ∇u 0 h 2 + ∇σ 0 h 2 dx ≤ C. (6.29)
In addition, we note the following result. Lemma 6.2. For p = 1, . . . , N P we have that

σ 0 min ξ 2 ≤ ξ T σ 0 h (P p ) ξ ≤ σ 0 max ξ 2 ∀ξ ∈ R d and tr(σ 0 h (P p )) ≤ b ⋆ . (6.30)
Proof. For the proof of the first result in (6.30), see the proof of Lemma 5.2 in [START_REF] Barrett | Existence and approximation of a (regularized) Oldroyd-B model[END_REF].

We now prove the second result in (6.30). On choosing χ = Iη, with η ∈ Q 1 h , in (6.28b) yields that z h := tr(σ

0 h ) -b ⋆ ∈ Q 1 h satisfies D [π h [z h η] + ∆t 0 ∇z h • ∇η] dx = D z η dx ∀η ∈ Q h 1 , (6.31)
where z := tr(σ 0 ) -b ⋆ ∈ L ∞ (D) and is non-positive on recalling (3.14)

. Choosing η = π h [z h ] + ∈ Q 1
h , it follows, on noting the Q 1 h version of (6.8) and (6.26) with g(

•) = [ • ] + , that D π h [z h ] + ] 2 + ∆t 0 ∇π h [z h ] + 2 dx ≤ d π h [z h ] 2 + + ∆t 0 ∇z h • ∇π h [z h ] + dx (6.32) = D z π h [z h ] + dx ≤ 0.
Hence π h [z h ] + ≡ 0 and so the second result in (6.30) holds. Furthermore, it follows from (6.28a,b), (6.29), (3.14), (6.3) and (6.9a,b) that, as h, ∆t 0 → 0 + ,

u 0 h → u 0 weakly in [L 2 (D)] d and σ 0 h → σ 0 weakly in [L 2 (D)] d×d . (6.33)
Our approximation (P ∆t α,δ,h ) of (P α,δ ) is then:

(P ∆t α,δ,h ) Setting (u 0 α,δ,h , σ 0 α,δ,h , ̺ 0 α,δ,h ) = (u 0 h , σ 0 h , tr(σ 0 h )) ∈ V 1 h × (S 1 h ∩ S >0,b ) × Q 1 h , with u 0 h and σ 0 h as defined in (6.28a,b), then for n = 1, . . . , N T find (u n α,δ,h , σ n α,δ,h , ̺ n α,δ,h ) ∈ V 1 h × S 1 h × Q 1 h such that for any test functions (v, φ, η) ∈ V 1 h × S 1 h × Q 1 h D Re u n α,δ,h -u n-1 α,δ,h ∆t n • v + Re 2 (u n-1 α,δ,h • ∇)u n α,δ,h • v -u n α,δ,h • (u n-1 α,δ,h • ∇)v (6.34a) + (1 -ε)∇u n α,δ,h : ∇v + ε Wi π h κ δ (σ n α,δ,h , ̺ n α,δ,h ) A δ (σ n α,δ,h , ̺ n α,δ,h ) β δ (σ n α,δ,h ) : ∇v dx = f n , v H 1 0 (D) , D π h σ n α,δ,h -σ n-1 α,δ,h ∆t n : φ + A δ (σ n α,δ,h , ̺ n α,δ,h ) β δ (σ n α,δ,h ) : φ Wi dx (6.34b) + D α∇σ n α,δ,h :: ∇φ -2∇u n α,δ,h : π h [κ δ (σ n α,δ,h , ̺ n α,δ,h ) φ β δ (σ n α,δ,h ) dx - D d m=1 d p=1 [u n-1 α,δ,h ] m Λ δ,m,p (σ n α,δ,h ) : ∂φ ∂x p dx = 0, D π h   ̺ n α,δ,h -̺ n-1 α,δ,h ∆t n η + tr A δ (σ n α,δ,h , ̺ n α,δ,h ) β δ (σ n α,δ,h ) Wi η   dx (6.34c) + D α∇̺ n α,δ,h • ∇η -2∇u n α,δ,h : π h [κ δ (σ n α,δ,h , ̺ n α,δ,h ) η β δ (σ n α,δ,h )] dx + b D d m=1 d p=1 [u n-1 α,δ,h ] m Λ δ,m,p 1 - ̺ n α,δ,h b ∂η ∂x p dx = 0.
In deriving (P ∆t α,δ,h ), we have noted (4.14), (1.5a), (6.19a,b) and (6.20). We note that on replacing A δ (σ n α,δ,h , tr(̺ n α,δ,h )) with I -G ′ δ (σ n α,δ,h ) and κ δ (σ n α,δ,h , ̺ n α,δ,h ) by 1 then (P ∆t α,δ,h ), (6.34a,b), collapses to the corresponding finite element approximation of Oldroyd-B with stress diffusion studied in [START_REF] Barrett | Existence and approximation of a (regularized) Oldroyd-B model[END_REF], see (5.34a,b) with no L cut-off there.

Before proving existence of a solution to (P ∆t α,δ,h ), we first derive a discrete analogue of the energy bound (5.8) for (P α,δ ).

Energy bound. On setting

F δ,h (v, φ, η) := Re 2 D v 2 dx - ε 2Wi D π h b G δ 1 - η b + tr (G δ (φ) + I) dx (6.35) ∀(v, φ, η) ∈ V 1 h × S 1 h × Q 1 h
, we have the following discrete analogue of Proposition 5.1.

Proposition 6.1. For n = 1, . . . , N T , a solution u n α,δ,h , σ n α,δ,h , ̺ n α,δ,h ∈ V 1 h ×S 1 h ×Q 1 h to (P ∆t α,δ,h ), (6.34a-c), if it exists, satisfies F δ,h (u n α,δ,h , σ n α,δ,h , ̺ n α,δ,h ) -F δ,h (u n-1 α,δ,h , σ n-1 α,δ,h , ̺ n-1 α,δ,h ) ∆t n (6.36) + Re 2∆t n D u n α,δ,h -u n-1 α,δ,h 2 dx + (1 -ε) D ∇u n α,δ,h 2 dx + ε 2Wi 2 D π h tr A δ (σ n α,δ,h , ̺ n α,δ,h ) 2 β δ (σ n α,δ,h ) dx + αεδ 2 2Wi D ∇π h [G ′ δ (σ n α,δ,h )] 2 + b ∇π h G ′ δ 1 - ̺ n α,δ,h b 2 dx ≤ f n , u n α,δ,h H 1 0 (D) ≤ 1 2 (1 -ε) D ∇u n α,δ,h 2 dx + 1 + C P 2(1 -ε) f n 2 H -1 (D) .
Proof. The proof is similar to that of Proposition 4.1, we choose as test functions

v = u n α,δ,h ∈ V 1 h , φ = -ε 2Wi π h [G ′ δ (σ n α,δ,h )] ∈ S 1 h and η = ε 2Wi π h G ′ δ 1 - ̺ n α,δ,h b ∈ Q 1 h in (6.34a-c),
f n , u n α,δ,h H 1 0 (D) ≥ F δ,h (u n α,δ,h , σ n α,δ,h , ̺ n α,δ,h ) -F δ,h (u n-1 α,δ,h , σ n-1 α,δ,h , ̺ n-1 α,δ,h ) ∆t n (6.37) + Re 2∆t n D u n α,δ,h -u n-1 α,δ,h 2 dx + (1 -ε) D ∇u n α,δ,h 2 dx + ε 2Wi 2 D π h tr A δ (σ n α,δ,h , ̺ n α,δ,h ) 2 β δ (σ n α,δ,h ) dx + αεδ 2 2Wi D ∇π h [G ′ δ (σ n α,δ,h )] 2 + b ∇π h G ′ δ 1 - ̺ n α,δ,h b 2 dx + ε 2Wi D u n-1 α,δ,h • ∇π h tr(H δ (G ′ δ (σ n α,δ,h ))) + bH δ G ′ δ 1 - ̺ n α,δ,h b dx.
The first desired inequality in (6.36) follows immediately from (6.37) on noting (6.1a,d), (1.12) and that π h : C(D) → Q 1 h . The second inequality in (6.36) follows immediately from (2.9) with

ν 2 = (1 -ε)/(1 + C P ). 6.4. Existence of discrete solutions. Proposition 6.2. Given (u n-1 α,δ,h , σ n-1 α,δ,h , ̺ n-1 α,δ,h ) ∈ V 1 h × S 1 h × Q 1 h such that D [tr(σ n-1 α,δ,h ) -̺ n-1 α,δ,h
] dx = 0 and for any time step ∆t n > 0, then there exists at least one solution

u n α,δ,h , σ n α,δ,h , ̺ n α,δ,h ∈ V 1 h × S 1 h × Q 1 h to (P ∆t α,δ,h ), (6.34a-c), such that D [tr(σ n α,δ,h ) -̺ n α,δ,h ] dx = 0. Proof.
The proof is similar to that of Proposition 4.2. We introduce the following inner product on the Hilbert space

V 1 h × S 1 h × Q 1 h ((w, ψ, ξ), (v, φ, η)) h D = D [w • v + π h [ψ : φ + ξ η] ] dx ∀(w, ψ, ξ), (v, φ, η) ∈ V 1 h × S 1 h × Q 1 h . Given (u n-1 α,δ,h , σ n-1 α,δ,h , ̺ n-1 α,δ,h ) ∈ V 1 h × S 1 h × Q 1 h , let F h : V 1 h × S 1 h × Q 1 h → V 1 h × S 1 h × Q 1 h be such that for any (w, ψ, ξ) ∈ V 1 h × S 1 h × Q 1 h F h (w, ψ, ξ), (v, φ, η) h D (6.38) := D Re w -u n-1 α,δ,h ∆t n • v + (1 -ε)∇w : ∇v + ε Wi π h [κ δ (ψ, ξ) A δ (ψ, ξ) β δ (ψ)] : ∇v + Re 2 (u n-1 α,δ,h • ∇)w • v -w • (u n-1 α,δ,h • ∇)v + α [∇ψ :: ∇φ + ∇ξ • ∇η] -2∇w : π h [ κ δ (ψ, ξ) [φ + ηI] β δ (ψ)] dx + D π h ψ -σ n-1 α,δ,h ∆t n : φ + A δ (ψ, ξ) β δ (ψ) Wi : φ + ξ -̺ n-1 α,δ,h ∆t n η + tr (A δ (ψ, ξ) β δ (ψ)) Wi η dx - D d m=1 d p=1 [u n-1 α,δ,h ] m Λ δ,m,p (ψ) : ∂φ ∂x p -b Λ δ,m,p 1 - ξ b ∂η ∂x p dx -f n , v H 1 0 (D) ∀(v, φ, η) ∈ V 1 h × S 1 h × Q 1 h . A solution (u n α,δ,h , σ n α,δ,h , ̺ n α,δ,h
) to (4.13a-c), if it exists, corresponds to a zero of F h . On recalling (6.19a,b), (6.18a,b) and (6.14), it is easily deduced that the mapping F h is continuous. For any

(w, ψ, ξ) ∈ V 1 h × S 1 h × Q 1 h , on choosing (v, φ, η) = w, -ε 2Wi π h [G ′ δ (ψ)], ε 2Wi π h G ′ δ 1 -ξ b
, we obtain analogously to (6.36) that

F h (w, ψ, ξ), w, - ε 2Wi π h [G ′ δ (ψ)], ε 2Wi π h G ′ δ 1 - ξ b h D (6.39) ≥ F δ,h (w, ψ, ξ) -F δ,h (u n-1 α,δ,h , σ n-1 α,δ,h , ̺ n-1 α,δ,h ) ∆t n + Re 2∆t n D w -u n-1 α,δ,h 2 dx + 1 -ε 2 D ∇w 2 dx - 1 + C P 2(1 -ε) f n 2 H -1 (D) + ε 2Wi 2 D π h tr (A δ (ψ, ξ)) 2 β δ (ψ) dx + αεδ 2 2Wi D ∇π h [G ′ δ (ψ)] 2 + b ∇π h G ′ δ 1 - ξ b 2 dx. Let (v, φ, η) h D := ((v, φ, η), (v, φ, η)) h D 1 2 = D v 2 + π h [ φ 2 + |η| 2 ] dx 1 2
.

If for any γ ∈ R >0 , the continuous mapping F h has no zero (u n α,δ,h , σ n α,δ,h , ̺ n α,δ,h ), which lies in the ball

B h γ := (v, φ, η) ∈ V 1 h × S 1 h × Q 1 h : (v, φ, η) h D ≤ γ ; then for such γ, we can define the continuous mapping G h γ : B h γ → B h γ such that for all (v, φ, η) ∈ B h γ G h γ (v, φ, η) := -γ F h (v, φ, η) F h (v, φ, η) h D
.

By the Brouwer fixed point theorem, G h γ has at least one fixed point (w γ , ψ γ , ξ γ ) in B h γ . Hence it satisfies (6.40) (w γ , ψ γ , ξ γ )

h D = G h γ (w γ , ψ γ , ξ γ ) h D = γ. In addition, (w γ , ψ γ , ξ γ ), (v, φ, η) h D = G h (w γ , ψ γ , ξ γ ), (v, φ, η) h D with (v, φ, η) = (0, I, -1) yields that D [tr(ψ γ ) -ξ γ ] dx = D [tr(σ n-1 α,δ,h ) -̺ n-1 α,δ,h ] dx = 0. (6.41)
On noting (6.4), we have that there exists a µ h ∈ R >0 such that for all φ ∈ S 1 h , (6.42)

π h [ φ ] 2 L ∞ (D) = π h [ φ 2 ] L ∞ (D) ≤ µ 2 h D π h [ φ 2 ] dx,
and an equivalent result holding for all η ∈ Q 1 h . It follows from (6.35), (6.41), (3.5), (3.7a), (6.42) and (6.40) that

F δ,h (w γ , ψ γ , ξ γ ) (6.43) = Re 2 D w γ 2 dx + ε 2Wi D π h tr(ψ γ -G δ (ψ γ ) -I) -b G δ 1 - ξ γ b -ξ γ dx ≥ Re 2 D w γ 2 dx + ε 4Wi D π h [ ψ γ + |ξ γ | ] dx -(2d + 3b)|D| ≥ Re 2 D w γ 2 dx - ε(2d + 3b)|D| 4Wi + ε 4Wi µ h γ π h [ ψ γ ] L ∞ (D) D π h [ ψ γ ] dx + π h [ |ξ γ | ] L ∞ (D) D π h [ |ξ γ | ] dx ≥ min Re 2 , ε 4Wi µ h γ D w γ 2 + π h [ ψ γ 2 + |ξ γ | 2 ] dx - ε(2d + 3b)|D| 4Wi = min Re 2 , ε 4Wi µ h γ γ 2 - ε(2d + 3b)|D| 4Wi .
Hence for all γ sufficiently large, it follows from (6.39) and (6.43) that (6.44)

F h (w γ , ψ γ , ξ γ ), w γ , - ε 2Wi π h [G ′ δ (ψ γ )], ε 2Wi π h G ′ δ 1 - ξ γ b h D ≥ 0.
On the other hand as (w γ , ψ γ , ξ γ ) is a fixed point of G h γ , we have that

F h (w γ , ψ γ , ξ γ ), w γ , - ε 2Wi π h [G ′ δ (ψ γ )], ε 2Wi π h G ′ δ 1 - ξ γ b h D (6.45) = - F h (w γ , ψ γ , ξ γ ) h D γ D w γ 2 - ε 2Wi π h ψ γ : G ′ δ (ψ γ ) -ξ γ G ′ δ 1 - ξ γ b dx.
It follows from (6.41), (3.5), (3.7b), and similarly to (6.43), on noting (6.42) and ( 6.40), that

D w γ 2 - ε 2Wi π h ψ γ : G ′ δ (ψ γ ) -ξ γ G ′ δ 1 - ξ γ b dx (6.46) = D w γ 2 + ε 2Wi π h ψ γ : (I -G ′ δ (ψ γ )) + ξ γ G ′ δ 1 - ξ γ b -1 dx ≥ D w γ 2 + ε 4Wi π h [ ψ γ + |ξ γ | ] -2(d + b) dx ≥ min 1, ε 4Wi µ h γ γ 2 - ε(d + b)|D| 2Wi .
Therefore on combining (6.45) and ( 6.46), we have for all γ sufficiently large that

F h (w γ , ψ γ , ξ γ ), w γ , - ε 2Wi π h [G ′ δ (ψ γ )], ε 2Wi π h G ′ δ 1 - ξ γ b h D < 0 , (6.47)
which obviously contradicts (6.44). Hence the mapping F h has a zero, (u n α,δ,h , σ n α,δ,h , ̺ n α,δ,h ) ∈ B h γ for γ sufficiently large. Finally, similarly to (6.41), it follows, on choosing (v, φ, η) = (0, I, -1)

in F h (u n α,δ,h , σ n α,δ,h , ̺ n α,δ,h ), (v, φ, η) h D = 0, that D [tr(σ n α,δ,h ) -̺ n α,δ,h ] dx = D [tr(σ n-1 α,δ,h ) - ̺ n-1 α,δ,h ] dx = 0.
We now have the analogue of Theorem 4.1.

Theorem 6.1. For any δ ∈ (0, min{ 1 2 , b}], N T ≥ 1 and any partitioning of [0, T ] into N T time steps, there exists a solution

{(u n α,δ,h , σ n α,δ,h , ̺ n α,δ,h )} NT n=1 ∈ [V 1 h ×S 1 h ×Q 1 h ] NT to (P ∆t α,δ,h
), (6.34a-c). In addition, it follows for n = 1, . . . , N T that D [tr(σ n α,δ,h ) -̺ n α,δ,h ] dx = 0 and

F δ,h (u n α,δ,h , σ n α,δ,h , ̺ n α,δ,h ) + 1 2 n m=1 D Re u m α,δ,h -u m-1 α,δ,h 2 + (1 -ε)∆t m ∇u m α,δ,h 2 dx (6.48) + ε 2Wi 2 n m=1 ∆t m D π h tr A δ (σ m α,δ,h , ̺ m α,δ,h ) 2 β δ (σ m α,δ,h ) dx + αεδ 2 2Wi n m=1 ∆t m D ∇π h [G ′ δ (σ m α,δ,h )] 2 + b ∇π h G ′ δ 1 - ̺ m α,δ,h b 2 dx ≤ F δ,h (u 0 h , σ 0 h , tr(σ 0 h )) + 1 + C P 2(1 -ε) n m=1 ∆t m f m 2 H -1 (D) ≤ C, which yields that max n=0,...,NT D u n α,δ,h 2 + σ n α,δ,h + |̺ n α,δ,h | + δ -1 π h [σ n α,δ,h ] -+ [b -̺ n α,δ,h ] - dx ≤ C; (6.49)
where C is independent of α, as well as δ, h and ∆t.

Proof. Existence and the stability result (6.48) follow immediately from Propositions 6.2 and 6.1, respectively, on noting (6.35), (3.1), (6.29), (6.30), (4.11a) and (3.14). The bounds (6.49) follow immediately from (6.48), on noting (3.6), (3.4b), (6.35), (3.5) and that

D π h b 1 - ̺ n α,δ,h b -G δ 1 - ̺ n α,δ,h b + tr(σ n α,δ,h -G δ (σ n α,δ,h )) dx (6.50) = D π h b 1 -G δ 1 - ̺ n α,δ,h b -tr(G δ (σ n α,δ,h )) dx ≤ C.
Remark 6.1. We recall that we have used S 1 h for the approximation of σ α,δ in (P ∆t α,δ,h ), (6.34a-c), due to the presence of the diffusion term in (5.5c). Secondly, due to the advective term in (5.5c), one has to introduce the variable ̺ n α,δ,h and its equation (6.34c) in (P ∆t α,δ,h ) in order to obtain the entropy bound (6.36). However, we now have a bound on π h [[b -̺ n α,δ,h ] -] in (6.49), as opposed to [b -tr(σ n δ,h )] -in (4.37). Now, it does not seem possible to pass to the limit δ → 0 in (P ∆t α,δ,h ) to prove well-posedness of the corresponding direct approximation of (P α ), i.e. (P ∆t α,h ) without the regularization δ, as we did for (P ∆t δ,h ) in subsection 4.5. Finally, we note the following Lemmas for later purposes.

Lemma 6.3. For all K k ∈ T h and φ ∈ [C(K k )] d×d S , we have, for r ∈ [1, ∞), that K k [π h [ φ r ] + |π h [ φ ]| r ] dx ≤ C K k π h φ r dx. (6.51) Proof. It follows immediately from (6.4) that K k [π h [ φ r ] + |π h [ φ ]| r ] dx ≤ 2 |K k | π h φ r L ∞ (K k ) ≤ C K k π h φ r dx,
and hence the desired result(6.51).

Lemma 6.4. Let g ∈ C 0,1 (R) with Lipschitz constant g Lip . For all K k ∈ T h , and for all q ∈ Q 1 h , φ ∈ S 1 h we have that

K k π h [β δ (φ)] -β δ (φ) 2 dx + max m,p=1,...,d K k Λ δ,m,p (φ) -β δ (φ) δ mp 2 dx (6.52a) ≤ C h 2 K k ∇φ 2 dx, K k π h [g(q)] -g(q) 2 dx ≤ C g 2 Lip h 2 K k ∇q 2 dx (6.52b) and K k π h [g(φ)] -g(φ) 2 dx ≤ C g 2 Lip h 2 K k ∇φ 2 dx.
In addition, if g is monotonic then, for all K k ∈ T h and for all q ∈ Q 1 h , we have that

K k π h [g(q)] -g(q) 2 dx ≤ C h 2 K k ∇π h [g(q)] 2 dx. (6.53)
Proof. The results (6.52a) are proved in Lemma 5.3 of [START_REF] Barrett | Existence and approximation of a (regularized) Oldroyd-B model[END_REF] for the case when β δ , Λ δ,m,p and S 1 h are replaced by β, Λ m,p and S 1 h ∩ S >0 . The proofs given there are trivially adapted to the present case. In fact, the proof of the first result in (6.52a) in [START_REF] Barrett | Existence and approximation of a (regularized) Oldroyd-B model[END_REF] is easily adapted to any function g ∈ C 0,1 (R). Hence, we have the results (6.52b).

The result (6.53) is a simple variation of (6.52b) and follows on noting that (6.54) where {P k j } d j=0 are the vertices of K k .

π h [g(q)] -g(q) L ∞ (K k ) ≤ max i, j=0,...,P k d |g(q(P k i )) -g(q(P k j ))|,
7. Convergence of (P ∆t α,δ,h ) to (P α ) in the case d = 2 7.1. Stability. Before proving our stability results, we introduce some further notation. We require the

L 2 projector R h : V → V 1 h defined by D (v -R h v) • w dx = 0 ∀w ∈ V 1 h . (7.1)
In addition, we require

P h : [L 2 (D)] d×d S → S 1 h defined by D π h [P h χ : φ] dx = D χ : φ dx ∀φ ∈ S 1 h . (7.2)
It is easily deduced for p = 1, . . . , N P and i, j = 1, . . . , d that

[P h χ] ij (P p ) = 1 D η p D [P h χ] ij η p dx, (7.3) where η p ∈ Q 1
h is such that η p (P r ) = δ pr for p, r = 1, . . . , N P . It follows from (7.2) and (6.7b) with φ = P h χ, in both cases, that

D P h χ 2 dx ≤ D π h [ P h χ 2 ] dx ≤ D χ 2 dx ∀χ ∈ [L 2 (D)] d×d S . (7.4)
We shall assume from now on that D is convex and that the family h } h>0 is quasi-uniform, i.e. h k ≥ C h, k = 1, . . . , N K . It then follows that

R h v H 1 (D) ≤ C v H 1 (D)
∀v ∈ V, (7.5) see Lemma 4.3 in Heywood and Rannacher [START_REF] Heywood | Finite element approximation of the nonstationary Navier-Stokes problem I: Regularity of solutions and second-order error estimates for spatial discretization[END_REF]. Similarly, it is easily established that

P h χ H 1 (D) ≤ C χ H 1 (D) ∀χ ∈ [H 1 (D)] d×d S . (7.6)
We also require the scalar analogue of P h , where d = 1 and S 1 h is replaced by Q 1 h , satisfying the analogues of (7.2)-(7.4) and (7.6).

Let

S

. We note that

χ, Eχ H 1 (D) = Eχ 2 H 1 (D) ∀χ ∈ ([H 1 (D)] d×d S ) ′ , (7.8) and E • H 1 (D) is a norm on ([H 1 (D)] d×d

S

) ′ . We also employ this operator in the scalar case, d = 1, i.e. E :

H 1 (D) ′ → H 1 (D).
Let V ′ be the topological dual of V with the space of weakly divergent free functions in [L 2 (D)] d being the pivot space. Let S : V ′ → V be such that Sw is the unique solution to the Helmholtz-Stokes problem

D [∇(Sw) : ∇v + (Sw) • v] dx = w, v V ∀v ∈ V, (7.9)
where •, • V denotes the duality pairing between V ′ and V. We note that

w, Sw V = Sw 2 H 1 (D) ∀w ∈ V ′ , (7.10) and S • H 1 (D) is a norm on the reflexive space V ′ .
We recall the following well-known Gagliardo-Nirenberg inequality.

Let r ∈ [2, ∞) if d = 2, and r ∈ [2, 6] if d = 3 and θ = d( 1 2 -1 r ). Then, there exists a positive constant C(D, r, d) such that η L r (D) ≤ C(D, r, d) η 1-θ L 2 (D) η θ H 1 (D)
∀η ∈ H 1 (D). (7.11) We recall also the following compactness results. Let Y 0 , Y and Y 1 be real Banach spaces, Y i , i = 0, 1, reflexive, with a compact embedding Y 0 ֒→ Y and a continuous embedding Y ֒→ Y 1 . Then, for µ i > 1, i = 0, 1, the following embedding is compact:

{ η ∈ L µ0 (0, T ; Y 0 ) : ∂η ∂t ∈ L µ1 (0, T ; Y 1 ) } ֒→ L µ0 (0, T ; Y); (7.12)
see Theorem 2.1 on p184 in Temam [START_REF] Temam | Navier-Stokes Equations. Theory and Numerical Analysis[END_REF]. Let X 0 , X and X 1 be real Hilbert spaces with a compact embedding X 0 ֒→ X and a continuous embedding X ֒→ X 1 . Then, for γ > 0, the following embedding is compact:

{ η ∈ L 2 (0, T ; X 0 ) : D γ t η ∈ L 2 (0, T ; X 1 ) } ֒→ L 2 (0, T ; X ), (7.13)
where D γ t η is the time derivative of order γ of η, which can be defined in terms of the Fourier transform of η; see Theorem 2.2 on p186 in Temam [START_REF] Temam | Navier-Stokes Equations. Theory and Numerical Analysis[END_REF].

Finally, we recall the discrete Gronwall inequality:

(r 0 ) 2 + (s 0 ) 2 ≤ (q 0 ) 2 , (7.14) (r m ) 2 + (s m ) 2 ≤ m-1 n=0 (η n ) 2 (r n ) 2 + m n=0 (q n ) 2 m ≥ 1 ⇒ (r m ) 2 + (s m ) 2 ≤ exp( m-1 n=0 (η n ) 2 ) m n=0 (q n ) 2 m ≥ 1.
Theorem 7.1. Under the assumptions of Theorem 6.1, there exists a solution

{(u n α,δ,h , σ n α,δ,h , ̺ n α,δ,h )} NT n=1 ∈ [V 1 h × S 1 h × Q 1 h ] NT of (P ∆t α,δ,h ), (6.34a-c)
, such that the bounds (6.48) and (6.49) hold.

Moreover, if d = 2, (6.27) holds and ∆t ≤ C ⋆ (ζ -1 ) α 1+ζ h 2 , for a ζ > 0 and a C ⋆ (ζ -1 ) ∈ R >0 sufficiently small, then the following bounds hold:

max n=0,...,NT D π h [ σ n α,δ,h 2 ] dx + NT n=1 D ∆t n α ∇σ n α,δ,h 2 + π h [ σ n α,δ,h -σ n-1 α,δ,h 2 ] dx ≤ C, (7.15a)
max n=0,...,NT D π h [ |̺ n α,δ,h | 2 ] dx + NT n=1 D ∆t n α ∇̺ n α,δ,h 2 + π h [ |̺ n α,δ,h -̺ n-1 α,δ,h | 2 ] dx ≤ C, (7.15b)
NT n=1 ∆t n 2 m=1 2 p=1 Λ δ,m,p (σ n α,δ,h ) 4 L 4 (D) + Λ δ,m,p 1 - ̺ n α,δ,h b 4 L 4 (D) dx ≤ C. (7.15c)
Proof. Existence and the bounds (6.48) and (6.49) were proved in Theorem 6.1.

On choosing φ = σ n α,δ,h in (6.34b), it follows from (3.2), and (4.15), on applying a Young's inequality, for any ζ > 0 that

1 2 D π h [ σ n α,δ,h 2 + σ n α,δ,h -σ n-1 α,δ,h 2 ] dx + ∆t n α D ∇σ n α,δ,h 2 dx (7.16)
+ ∆t n Wi D π h G ′ δ 1 - ̺ n α,δ,h b tr (β δ (σ n α,δ,h )) 2 dx ≤ 1 2 D π h [ σ n-1 α,δ,h 2 ] dx + 2∆t n D ∇u n α,δ,h : π h [κ δ (σ n α,δ,h , ̺ n α,δ,h ) σ n α,δ,h β δ (σ n α,δ,h )] dx + ∆t n Wi D π h G ′ δ 1 - ̺ n α,δ,h b tr β δ (σ n α,δ,h ) (β δ (σ n α,δ,h ) -σ n α,δ,h ) dx + ∆t n Wi σ n α,δ,h L 1 (D) + ∆t n D d m=1 d p=1 [u n-1 α,δ,h ] m Λ δ,m,p (σ n α,δ,h ) : ∂σ n α,δ,h ∂x p dx ≤ 1 2 D π h [ σ n-1 α,δ,h 2 ] dx + ∆t n Wi σ n α,δ,h L 1 (D) + C δ -1 ∆t n D π h tr β δ (σ n α,δ,h ) β δ (σ n α,δ,h ) -σ n α,δ,h dx + 2 ∆t n ∇u n α,δ,h L 2 (D) π h [κ δ (σ n α,δ,h , ̺ n α,δ,h ) σ n α,δ,h β δ (σ n α,δ,h )] L 2 (D) + C ∆t n u n-1 α,δ,h L 2(2+ζ) ζ (D) Λ δ,m,p (σ n α,δ,h ) L 2+ζ (D) ∇σ n α,δ,h L 2 (D) .
We deduce from (3.3), (1.8), (6.7b), (1.5b), (5.2), (5.3), (6.51) and (7.11), as d = 2, that

tr β δ (σ n α,δ,h ) β δ (σ n α,δ,h ) -σ n α,δ,h ≤ δ tr β δ (σ n α,δ,h ) -σ n α,δ,h ≤ C δ (δ + σ n α,δ,h ), (7.17a) π h [κ δ (σ n α,δ,h , ̺ n α,δ,h ) σ n α,δ,h β δ (σ n α,δ,h )] 2 L 2 (D) (7.17b)
≤ D π h κ δ (σ n α,δ,h , ̺ n α,δ,h ) σ n α,δ,h β δ (σ n α,δ,h ) 2 dx ≤ C b D π h σ n α,δ,h 2 β δ (σ n α,δ,h ) dx ≤ C δ 3 + D π h [ σ n α,δ,h 3 ] dx ≤ C δ 3 + σ n α,δ,h 3
L 3 (D) ≤ C δ 3 + σ n α,δ,h 2
L 2 (D) σ n α,δ,h H 1 (D)
. Similarly, as d = 2, it follows from (6.24b), (3.3), (6.4) and (7.11)

that for all ζ > 0 Λ δ,m,p (σ n α,δ,h ) 2+ζ L 2+ζ (D) ≤ NK k=1 |K k | Λ δ,m,p (σ n α,δ,h ) 2+ζ L ∞ (K k ) ≤ C NK k=1 |K k | π h [β δ (σ n α,δ,h )] 2+ζ L ∞ (K k ) (7.18)
≤ C δ 2+ζ + σ n α,δ,h 2+ζ L 2+ζ (D) ≤ C + C(ζ) σ n α,δ,h 2
L 2 (D) σ n α,δ,h ζ H 1 (D)
. In addition, as d = 2, we note from (7.11), (1.13) and (6.49) that for all ζ > 0

u n-1 α,δ,h L 2(2+ζ) ζ (D) ≤ C(ζ -1 ) u n-1 α,δ,h ζ 2+ζ L 2 (D) u n-1 α,δ,h 2 2+ζ H 1 (D) ≤ C(ζ -1 ) ∇u n-1 α,δ,h
D π h [ σ n α,δ,h 2 + σ n α,δ,h -σ n-1 α,δ,h 2 ] dx + ∆t n α D ∇σ n α,δ,h 2 dx (7.20)
≤ D π h [ σ n-1 α,δ,h 2 ] dx + C(ζ -1 ) ∆t n α -(1+ζ) 1 + ∇u n α,δ,h 2
L 2 (D) + ∇u n-1 α,δ,h 2 L 2 (D) 1 + σ n α,δ,h 2
L 2 (D)
. Hence, summing (7.20) from n = 1, . . . , m for m = 1, . . . , N T yields, on noting (6.7b), that for any

ζ > 0 D π h [ σ m α,δ,h 2 ] dx + α m n=1 ∆t n D ∇σ n α,δ,h 2 dx + m n=1 D π h [ σ n α,δ,h -σ n-1 α,δ,h 2 ] dx (7.21) ≤ D π h [ σ 0 h 2 ] dx + C(ζ -1 ) α -(1+ζ) + C(ζ -1 ) α -(1+ζ) m n=1 ∆t n 1 + n k=n-1 ∇u k α,δ,h 2 L 2 (D) D π h [ σ n α,δ,h 2 ] dx.
Applying the discrete Gronwall inequality (7.14) to (7.21), and noting (6.27), (6.8), (6.29), (6.5b), (6.48), (6.49), and that ∆t

≤ C ⋆ (ζ -1 ) α 1+ζ h 2 , for a ζ > 0 where C ⋆ (ζ -1
) is sufficiently small, yields the bounds (7.15a).

Similarly to (7.16) on choosing η = ̺ n α,δ,h in (6.34c), it follows from (4.15), (1.3b), (1.4b) and (1.8), on applying a Young's inequality, for any ζ > 0 that

1 2 D π h [ |̺ n α,δ,h | 2 + |̺ n α,δ,h -̺ n-1 α,δ,h | 2 ] dx + ∆t n α D ∇̺ n α,δ,h 2 dx (7.22)
≤ 1 2 D π h [ |̺ n-1 α,δ,h | 2 ] dx + 2∆t n D ∇u n α,δ,h : π h [κ δ (σ n α,δ,h , ̺ n α,δ,h ) ̺ n α,δ,h β δ (σ n α,δ,h )] dx - ∆t n Wi D π h tr A δ (σ n α,δ,h , ̺ n α,δ,h ) β δ (σ n α,δ,h ) ̺ n α,δ,h dx -∆t n b D d m=1 d p=1 [u n-1 α,δ,h ] m Λ δ,m,p 1 - ̺ n α,δ,h b ∂̺ n α,δ,h ∂x p dx ≤ 1 2 D π h [ |̺ n-1 α,δ,h | 2 ] dx + 2 ∆t n ∇u n α,δ,h L 2 (D) π h [κ δ (σ n α,δ,h , ̺ n α,δ,h ) ̺ n α,δ,h β δ (σ n α,δ,h )] L 2 (D) + C ∆t n u n-1 α,δ,h L 2(2+ζ) ζ (D) Λ δ,m,p 1 - ̺ n α,δ,h b L 2+ζ (D) ∇̺ n α,δ,h L 2 (D) + ∆t n 2Wi D π h tr A δ (σ n α,δ,h , ̺ n α,δ,h ) 2 β δ (σ n α,δ,h ) + d 1 2 (̺ n α,δ,h ) 2 β δ (σ n α,δ,h ) dx.
Similarly to (7.17b), as d = 2, we deduce from (6.7b), (5.2), (5.3), (3.3), (6.4), (7.11) and (7.15a) that

π h [κ δ (σ n α,δ,h , ̺ n α,δ,h ) ̺ n α,δ,h β δ (σ n α,δ,h )] 2 L 2 (D) ≤ C b D π h (̺ n α,δ,h ) 2 β δ (σ n α,δ ,h ) dx (7.23)
≤ C 1 + σ n α,δ,h L 3 (D) ̺ n α,δ,h 2 L 3 (D) ≤ C 1 + σ n α,δ,h 1 3 H 1 (D) ̺ n α,δ,h 4 3 L 2 (D) ̺ n α,δ,h 2 3
H 1 (D) . Similarly to (7.18), as d = 2, it follows from (6.24a), (3.3), (6.4) and (7.11)

that for all ζ > 0 Λ δ,m,p 1 - ̺ n α,δ,h b 2+ζ L 2+ζ (D) ≤ C + C(ζ) ̺ n α,δ,h 2
L 2 (D) ̺ n α,δ,h ζ H 1 (D) . ( 7
D π h [ |̺ n α,δ,h 2 + |̺ n α,δ,h -̺ n-1 α,δ,h | 2 ] dx + ∆t n α D ∇̺ n α,δ,h 2 dx (7.25) ≤ D π h [ |̺ n-1 α,δ,h | 2 ] dx + 2 ∆t n ∇u n α,δ,h 2 L 2 (D) + ∆t n Wi D π h tr A δ (σ n α,δ,h , ̺ n α,δ,h ) 2 β δ (σ n α,δ,h ) dx + C(ζ -1 ) ∆t n α -(1+ζ) 1 + σ n α,δ,h 2
H 1 (D) + ∇u n-1 α,δ,h 2 L 2 (D) 1 + ̺ n α,δ,h 2
L 2 (D)
. Hence, summing (7.25) from n = 1, . . . , m for m = 1, . . . , N T yields, on noting (6.7b) and (6.48), that for any ζ > 0

D π h [ |̺ m α,δ,h | 2 ] dx + α m n=1 ∆t n D ∇̺ n α,δ,h 2 dx + m n=1 D π h [ |̺ n α,δ,h -̺ n-1 α,δ,h 2 ] dx (7.26) ≤ D π h [ (tr(σ 0 h )) 2 ] dx + C(ζ -1 ) α -(1+ζ) + C(ζ -1 ) α -(1+ζ) m n=1 ∆t n 1 + σ n α,δ,h 2
H 1 (D) + ∇u n-1 α,δ,h 2 L 2 (D) D π h [ |̺ n α,δ,h | 2 ] dx.
Applying the discrete Gronwall inequality (7.14) to (7.26), and noting (6.27), (6.30), (6.5a,b), (6.48), (6.49), (6.8), (7.15a), and that ∆t ≤ C ⋆ (ζ -1 ) α 1+ζ h 2 , for a ζ > 0 where C ⋆ (ζ -1 ) is sufficiently small, yields the bounds (7.15b). Finally, the desired result (7.15c) follows immediately from (7.18) and (7.24) with ζ = 2, on noting (7.15a,b) and (6.8).

Remark 7.1. Our final convergence result will be restricted to d = 2 for the same reason as why our result for Oldroyd-B in [START_REF] Barrett | Existence and approximation of a (regularized) Oldroyd-B model[END_REF] was restricted to d = 2. For example, the control of the term (7.17b) necessitates the restriction to d = 2. This could be overcome for FENE-P by replacing the regularization β δ in (P α,δ ), (5.5a-h), by β b δ and using test functions based on G b δ in place of G δ , where (G b δ ) ′ (s) = β b δ (s) for all s ∈ R. On making a similar change to our numerical approximation (P h,∆t α,δ ), (6.34a-c), it is then possible to prove the analogues of Theorem 6.1 and Theorem 7.1 for d = 2 and 3 with no restriction on ∆t. One can then establish analogues of Lemmas 7.1, 7.2 and Theorem 7.2, below, but now the limits involve the cut-off b, i.e. u b α , σ b α and ̺ b α . Moreover, it does not seem possible to establish that ̺ b α = tr(σ b α ), as we now have Λ b δ,m,p in place of Λ δ,m,p in the analogue of (7.60). Hence, in taking the limit δ, h, ∆t → 0 + , the last term in the analogue of (7.61) would have b β b 1 -

̺ b α b + tr(β b (σ b α )) instead of -(̺ b α -tr(σ b α ))
, where β b (s) = min{s, b}. Therefore, the extension of the convergence analysis in this paper to d = 3 and with a weaker restriction on ∆t will be a topic of further research. Lemma 7.1. Under all of the assumptions of Theorem 7.1, the solution {(u n α,δ,h , σ n α,δ,h , ̺ n α,δ,h )} NT n=1 of (P ∆t α,δ,h ), (6.34a-c), satisfies the following bounds:

NT n=1 ∆t n S u n α,δ,h -u n-1 α,δ,h ∆t n 4 ϑ H 1 (D) ≤ C, (7.27a) NT n=1 ∆t n E σ n α,δ,h -σ n-1 α,δ,h ∆t n 2 H 1 (D) + NT n=1 ∆t n E ̺ n α,δ,h -̺ n-1 α,δ,h ∆t n 2 H 1 (D) ≤ C, (7.27b) NT n=1 ∆t n π h κ δ (σ n α,δ,h , ̺ n α,δ,h ) A δ (σ n α,δ,h , ̺ n α,δ,h ) β δ (σ n α,δ,h ) 2 L 2 (D) ≤ C, (7.27c) NT n=1 ∆t n π h A δ (σ n α,δ,h , ̺ n α,δ,h ) β δ (σ n α,δ,h ) 8 5 L 8 5 (D) ≤ C, (7.27d)
where ϑ ∈ [START_REF] Barrett | Convergence of a finite-element approximation of surfactant spreading on a thin film in the presence of van der Waals forces[END_REF][START_REF] Barrett | Existence and equilibration of global weak solutions to kinetic models for dilute polymers I: Finitely extensible nonlinear bead-spring chains[END_REF] and C in (7.27a,c) is independent of α, as well as δ, h and ∆t.

Proof. On choosing

v = R h S u n α,δ,h -u n-1 α,δ,h ∆tn ∈ V 1
h in (6.34a) yields, on noting (7.1), (7.10), (7.5) and Sobolev embedding, that

Re S u n α,δ,h -u n-1 α,δ,h ∆t n 2 H 1 (D) = Re D u n α,δ,h -u n-1 α,δ,h ∆t n • R h S u n α,δ,h -u n-1 α,δ,h ∆t n dx (7.28) = ε Wi D π h κ δ (σ n α,δ,h , ̺ n α,δ,h ) A δ (σ n α,δ,h , ̺ n α,δ,h ) β δ (σ n α,δ,h ) : ∇ R h S u n α,δ,h -u n-1 α,δ,h ∆t n dx -(1 -ε) D ∇u n α,δ,h : ∇ R h S u n α,δ,h -u n-1 α,δ,h ∆t n dx - Re 2 D (u n-1 α,δ,h • ∇)u n α,δ,h • R h S u n α,δ,h -u n-1 α,δ,h ∆t n dx + Re 2 D u n α,δ,h • (u n-1 α,δ,h • ∇) R h S u n α,δ,h -u n-1 α,δ,h ∆t n dx + f n , R h S u n α,δ,h -u n-1 α,δ,h ∆t n H 1 0 (D) ≤ C π h κ δ (σ n α,δ,h , ̺ n α,δ,h ) A δ (σ n α,δ,h , ̺ n α,δ,h ) β δ (σ n α,δ,h ) 2 L 2 (D) + ∇u n α,δ,h 2 L 2 (D) + u n-1 α,δ,h u n α,δ,h 2
L 2 (D) + u n-1 α,δ,h ∇u n α,δ,h 2
L 1+θ (D) + f n 2 H -1 (D) ,
where θ > 0 as d = 2. It follows from (6.7b) and (5.4) that

π h κ δ (σ n α,δ,h , ̺ n α,δ,h ) A δ (σ n α,δ,h , ̺ n α,δ,h ) β δ (σ n α,δ,h ) 2 L 2 (D) (7.29) ≤ D π h κ δ (σ n α,δ,h , ̺ n α,δ,h ) A δ (σ n α,δ,h , ̺ n α,δ,h ) β δ (σ n α,δ,h ) 2 dx ≤ b D π h tr A δ (σ n α,δ,h , ̺ n α,δ,h ) 2 β δ (σ n α,δ,h ) dx.
Applying the Cauchy-Schwarz and the algebraic-geometric mean inequalities, in conjunction with (7.11), for d = 2, and the Poincaré inequality (1.13) yields that

u n-1 α,δ,h u n α,δ,h 2
L 2 (D) ≤ u n-1 α,δ,h 2
L 4 (D) u n α,δ,h 2
L 4 (D) ≤ 1 2 n m=n-1 u m α,δ,h 4 L 4 (D) (7.30) ≤ C n m=n-1 u m α,δ,h 2
L 2 (D) ∇u m α,δ,h 2
L 2 (D) .
Similarly, we have for any θ ∈ (0, 1), as d = 2, but now using a Young's inequality

u n-1 α,δ,h ∇u n α,δ,h 2
L 1+θ (D) ≤ u n-1 α,δ,h 2 L 2(1+θ) 1-θ (D) ∇u n α,δ,h 2 L 2 (D) (7.31) ≤ C u n-1 α,δ,h 2(1-θ) 1+θ L 2 (D) ∇u n-1 α,δ,h 4θ 1+θ L 2 (D) ∇u n α,δ,h 2 L 2 (D) ≤ C u n-1 α,δ,h 2(1-θ) 1+θ L 2 (D) n m=n-1 ∇u m α,δ,h 2(1+3θ) 1+θ L 2 (D) .
On taking the 2 ϑ power of both sides of (7.28), multiplying by ∆t n , summing from n = 1, . . . , N T and noting (7.29), (7.30), (7.31) with θ = ϑ-2 6-ϑ ⇔ ϑ = 2(1+3θ) (1+θ) ∈ (2, 4), (6.27), (4.11a), (6.48), (6.49) and (6.29) yields that

NT n=1 ∆t n S u n α,δ,h -u n-1 α,δ,h ∆t n 4 ϑ H 1 (D) (7.32) ≤ C NT n=1 ∆t n D π h tr A δ (σ n α,δ,h , ̺ n α,δ,h ) 2 β δ (σ n α,δ,h ) dx 2 ϑ + C NT n=1 ∆t n ∇u n α,δ,h 2
L 2 (D) + f n 2 H -1 (D) 2 ϑ + C 1 + max n=0,...,NT u n α,δ,h 2 L 2 (D) NT n=0 ∆t n ∇u n α,δ,h 2 L 2 (D)
≤ C, and hence the bound (7.27a) with C independent of α.

Choosing φ = P h E σ n α,δ,h -σ n-1 α,δ,h ∆tn ∈ S 1
h in (6.34b) yields, on noting (7.2) and (7.8), that

E σ n α,δ,h -σ n-1 α,δ,h ∆t n 2 H 1 (D) (7.33) = D π h σ n α,δ,h -σ n-1 α,δ,h ∆t n : P h E σ n α,δ,h -σ n-1 α,δ,h ∆t n dx = - 1 Wi D π h A δ (σ n α,δ,h , ̺ n α,δ,h ) β δ (σ n α,δ,h ) : P h E σ n α,δ,h -σ n-1 α,δ,h ∆t n dx -α D ∇σ n α,δ,h :: ∇ P h E σ n α,δ,h -σ n-1 α,δ,h ∆t n dx + 2 D ∇u n α,δ,h : π h κ δ (σ n α,δ,h , ̺ n α,δ,h ) P h E σ n α,δ,h -σ n-1 α,δ,h ∆t n β δ (σ n α,δ,h ) dx + D d m=1 d p=1 [u n-1 α,δ,h ] m Λ δ,m,p (σ n α,δ,h ) : ∂ ∂x p P h E σ n α,δ,h -σ n-1 α,δ,h ∆t n dx.
It follows from (1.5a,b), (1.3b), (1.4b) and (1.8) that for any

ζ ∈ R >0 D π h A δ (σ n α,δ,h , ̺ n α,δ,h ) β δ (σ n α,δ,h ) : P h E σ n α,δ,h -σ n-1 α,δ,h ∆t n dx (7.34) ≤ D π h A δ (σ n α,δ,h , ̺ n α,δ,h ) [β δ (σ n α,δ,h )] 1 2 P h E σ n α,δ,h -σ n-1 α,δ,h ∆t n [β δ (σ n α,δ,h )] 1 2 dx ≤ ζ -1 D π h tr (A δ (σ n α,δ,h , ̺ n α,δ,h )) 2 β δ (σ n α,δ,h ) dx + ζ d 1 2 D π h   P h E σ n α,δ,h -σ n-1 α,δ,h ∆t n 2 β δ (σ n α,δ,h )   dx.
Similarly to (7.17b), it follows from (6.7b), (1.5b), (5.2), (5.3) and (1.8) that

π h κ δ (σ n α,δ,h , ̺ n α,δ,h ) P h E σ n α,δ,h -σ n-1 α,δ,h ∆t n β δ (σ n α,δ,h ) 2 L 2 (D) (7.35) ≤ C D π h   P h E σ n α,δ,h -σ n-1 α,δ,h ∆t n 2 β δ (σ n α,δ,h )   dx.
In addition, (3.3), (6.4), (6.8) and (7.15a) imply that for all φ ∈ S 1

h D π h φ 2 β δ (σ n α,δ,h ) dx ≤ NK k=1 σ n α,δ,h L ∞ (K k ) + δ K k φ 2 dx (7.36) ≤ C σ n α,δ,h L 2 (D) + δ φ 2 L 4 (D) ≤ C φ 2 H 1 (D)
. Combining (7.33)-(7.36), yields, on noting (7.6), that

E σ n α,δ,h -σ n-1 α,δ,h ∆t n 2 H 1 (D) (7.37) ≤ C D π h tr (A δ (σ n α,δ,h , ̺ n α,δ,h )) 2 β δ (σ n α,δ,h ) dx + α ∇σ n α,δ,h 2 L 2 (D) + ∇u n α,δ,h 2
L 2 (D) + u n-1 α,δ,h 2 L 4 (D) D 2 m=1 2 p=1 Λ δ,m,p (σ n α,δ,h ) 2 L 4 (D) dx .
Similarly to (7.33)-(7.37), choosing

η = P h E ̺ n α,δ,h -̺ n-1 α,δ,h ∆tn ∈ Q 1 h in (6.34c) yields that E ̺ n α,δ,h -̺ n-1 α,δ,h ∆t n 2 H 1 (D) (7.38) ≤ C D π h tr (A δ (σ n α,δ,h , ̺ n α,δ,h )) 2 β δ (σ n α,δ,h ) dx + α ∇̺ n α,δ,h 2 L 2 (D) + ∇u n α,δ,h 2
L 2 (D) + u n-1 α,δ,h 2 L 4 (D) D 2 m=1 2 p=1 Λ δ,m,p 1 - ̺ n α,δ,h b 2 L 4 (D)
dx .

Multiplying (7.37) and (7.38) by ∆t n , summing from n = 1, ..., N T and noting (6.48), (6.49), (7.30), (6.27), (6.29) and (7.15a-c) yields the bounds (7.27b). Multiplying (7.29) by ∆t n , summing from n = 1, ..., N T and noting (6.48) yields the result (7.27c) with C independent of α. Finally, it follows from (6.7a) that

π h A δ (σ n α,δ,h , ̺ n α,δ,h ) β δ (σ n α,δ,h ) 8 5 L 8 5 (D) (7.39)
≤ D π h A δ (σ n α,δ,h , ̺ n α,δ,h ) [β δ (σ n α,δ,h )] 1 2 2 4 5 π h [β δ (σ n α,δ,h )] 1 2 2 4 5 dx ≤ π h A δ (σ n α,δ,h , ̺ n α,δ,h ) [β δ (σ n α,δ,h )] 1 2 2 4 5 L 1 (D) π h [β δ (σ n α,δ,h )]
≤ NT n=1 ∆t n π h A δ (σ n α,δ,h , ̺ n α,δ,h ) [β δ (σ n α,δ,h )] 1 2 2 L 1 (D) 4 5 × NT n=1 ∆t n π h β δ (σ n α,δ,h ) 4 L 4 (D) 1 5 ≤ C NT n=1 ∆t n π h β δ (σ n α,δ,h ) 4 L 4 (D) 1 5 ≤ C δ 4 + NT n=1 ∆t n σ n α,δ,h 4
L 4 (D) 1 5 ≤ C δ 4 + NT n=1 ∆t n σ n α,δ,h 2
L 2 (D) σ n α,δ,h 2
H 1 (D) 1 5
≤ C.

Hence, we have the desired result (7.27d).

Unfortunately, the bound (7.27a) is not useful for obtaining compactness via (7.12), see the discussion in the proof of Theorem 7.2 below. Instead one has to exploit the compactness result (7.13). This we now do, by following the proof of Lemma 5.6 on p237 in Temam [START_REF] Temam | Navier-Stokes Equations. Theory and Numerical Analysis[END_REF]. Here the introduction of κ δ (σ n α,δ,h , ̺ n α,δ,h ) in the extra stress term, as discussed above (5.2), is crucial, as this yields an L 2 (D T ) bound on this stress term via the bound (5.4). For this purpose, we introduce the following notation in line with (4.10) as a collective symbol for u ∆t α,δ,h , u ∆t,± α,δ,h . We also define σ ∆t(,±) α,δ,h and ̺ ∆t(,±) α,δ,h similarly to (7.41a,b). Using the notation (7.41a,b), (P ∆t α,δ,h ), i.e. (6.34a-c) multiplied by ∆t n and summed for n = 1, . . . , N T , can be restated as:

. Let u ∆t α,δ,h ∈ C([0, T ]; V 1 h ) and u ∆t,± α,δ,h ∈ L ∞ (0, T ; V 1 h ) be such that for n = 1, . . . , N T u ∆t α,δ,h (t, •) := t -t n-1 ∆t n u n α,δ,h (•) + t n -t ∆t n u n-1 α,δ,h (•) t ∈ [t n-1 , t n ], (7.41a) u ∆t,+ α,δ,h (t, •) := u n α,δ,h (•), u ∆t,- α,δ,h (t, •) := u n-1 α,δ,h (•) t ∈ [t n-
DT Re ∂u ∆t α,δ,h ∂t • v + (1 -ε)∇u ∆t,+ α,δ,h : ∇v dx dt (7.43a) + Re 2 DT (u ∆t,- α,δ,h • ∇)u ∆t,+ α,δ,h • v -(u ∆t,- α,δ,h • ∇)v • u ∆t,+ α,δ,h dx dt + ε Wi DT π h κ δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) A δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) β δ (σ ∆t,+ α,δ,h ) : ∇v dx dt = T 0 f + , v H 1 0 (D) dt ∀v ∈ L 2 (0, T ; V
∂φ ∂x p dx dt = 0 ∀φ ∈ L 2 (0, T ; S 1 h ), DT π h   ∂̺ ∆t α,δ,h ∂t η + tr A δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) β δ (σ ∆t,+ α,δ,h ) Wi η   dx dt (7.43c) + α DT ∇̺ ∆t,+ α,δ,h • ∇η dx dt -2 DT ∇u ∆t,+ α,δ,h : π h [κ δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) η β δ (σ ∆t,+ α,δ,h )] dx dt + b DT d m=1 d p=1 [u ∆t,- α,δ,h ] m Λ δ,m,p 1 - ̺ ∆t,+ α,δ,h b ∂η ∂x p dx dt = 0 ∀η ∈ L 2 (0, T ; Q 1 h )
subject to the initial conditions u ∆t α,δ,h (0) = u 0 h , σ ∆t α,δ,h (0) = σ 0 h and ̺ ∆t α,δ,h (0) = tr(σ 0 h ).

Lemma 7.2. Under all of the assumptions of Theorem 7.1, the solution (u ∆t α,δ,h , σ ∆t α,δ,h , ̺ ∆t α,δ,h ) of (P ∆t α,δ,h ), (7.43a-c), satisfies the following bound:

T 0 D γ t u ∆t α,δ,h 2
L 2 (D) dt ≤ C, (7.44)
where γ ∈ (0, 1 4 ) and C is independent of α, as well as δ, h and ∆t.

Proof. Equation (7.43a) can be reinterpreted as

Re d dt D u ∆t α,δ,h (t) • v dx = D ∇g ∆t,+ (t) : ∇v dx ∀v ∈ V 1 h , t ∈ (0, T ), (7.45) where g ∆t,+ (t) ∈ V 1 h is defined by D ∇g ∆t,+ (t) : ∇v dx = f + (t), v H 1 0 (D) -(1 -ε) D ∇u ∆t,+ α,δ,h (t) : ∇v dx (7.46) - Re 2 D (u ∆t,- α,δ,h (t) • ∇)u ∆t,+ α,δ,h (t) • v -(u ∆t,- α,δ,h (t) • ∇)v • u ∆t,+ α,δ,h (t) dx - ε Wi D π h κ δ (σ ∆t,+ α,δ,h (t), ̺ ∆t,+ α,δ,h (t)) A δ (σ ∆t,+ α,δ,h (t), ̺ ∆t,+ α,δ,h (t)) β δ (σ ∆t,+ α,δ,h (t)) : ∇v dx.
Similarly to (7.28) and (7.29) with θ ∈ (0, 1), it follows from (7.46) that

∇g ∆t,+ (t) L 2 (D) ≤ C f + (t) H -1 (D) + ∇u ∆t,+ α,δ,h (t) L 2 (D) (7.47) + u ∆t,- α,δ,h (t) u ∆t,+ α,δ,h (t) L 2 (D) + u ∆t,- α,δ,h (t) ∇u ∆t,+ α,δ,h (t) L 1+θ (D) + D tr A δ (σ ∆t,+ α,δ,h (t), ̺ ∆t,+ α,δ,h (t)) 2 β δ (σ ∆t,+ α,δ,h (t)) dx 1 2
.

On noting (7.30), (7.31) and the bound on u n α,δ,h in (6.49), we deduce from (7.47) that

∇g ∆t,+ (t) L 2 (D) ≤ C 1 + f + (t) H -1 (D) + ∇u ∆t,- α,δ,h (t) 1+3θ 1+θ L 2 (D) + ∇u ∆t,+ α,δ,h (t) 1+3θ 1+θ L 2 (D) (7.48)
+ D tr A δ (σ ∆t,+ α,δ,h (t), ̺ ∆t,+ α,δ,h (t)) 2 β δ (σ ∆t,+ α,δ,h (t)) dx 1 2
.

Similarly to (7.32), on recalling (6.29), (4.11a) and (6.49), we deduce from (7.48) that T 0 ∇g ∆t,+ (t) where ϑ ∈ (2, 4] and C is independent of α, as well as δ, h and ∆t. The rest of the proof follows as on p1825-6 in [START_REF] Barrett | Existence and approximation of a (regularized) Oldroyd-B model[END_REF], which is based on the proof of Lemma 5.6 on p237 in [START_REF] Temam | Navier-Stokes Equations. Theory and Numerical Analysis[END_REF]. 7.2. Convergence. It follows from (6.48), (6.49), (7.15a-c), (6.29), (6.8), (7.27a-d), (7.44) and (7.41a-c) that sup t∈(0,T )

u ∆t(,±) α,δ,h 2 L 2 (D) + T 0 ∇u ∆t(,±) α,δ,h 2 L 2 (D) + u ∆t,+ α,δ,h -u ∆t,- α,δ,h 2 L 2 (D) ∆(t) dt ≤ C, (7.50a) sup t∈(0,T ) σ ∆t(,±) α,δ,h 2 L 2 (D) + T 0 α ∇σ ∆t(,±) α,δ,h 2 L 2 (D) + σ ∆t,+ α,δ,h -σ ∆t,- α,δ,h 2 L 2 (D) ∆(t) dt ≤ C, (7.50b) sup t∈(0,T ) ̺ ∆t(,±) α,δ,h 2 L 2 (D) + T 0 α ∇̺ ∆t(,±) α,δ,h 2 L 2 (D) + ̺ ∆t,+ α,δ,h -̺ ∆t,- α,δ,h 2 L 2 (D) ∆(t) dt (7.50c) + δ 2 T 0 ∇π h G ′ δ 1 - ̺ ∆t,+ α,δ,h b 2 L 2 (D) dt ≤ C, DT π h tr A δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) 2 β δ (σ ∆t,+ α,δ,h ) dx dt ≤ C, (7.50d) sup t∈(0,T ) D π h [σ ∆t(,±) α,δ,h ] -+ |[b -̺ ∆t(,±) α,δ,h ] -| dx ≤ C δ, (7.50e) T 0   S ∂u ∆t α,δ,h ∂t 4 ϑ H 1 (D) + D γ t u ∆t α,δ,h 2 L 2 (D)   dt ≤ C, (7.50f) T 0   E ∂σ ∆t α,δ,h ∂t 2 H 1 (D) + E ∂̺ ∆t α,δ,h ∂t 2 H 1 (D)   dt ≤ C, (7.50g) π h κ δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) A δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) β δ (σ ∆t,+ α,δ,h ) L 2 (DT ) ≤ C, (7.50h) π h A δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) β δ (σ ∆t,+ α,δ,h ) L 8 5 (DT ) ≤ C, (7.50i)
where ϑ ∈ (2, 4], γ ∈ (0, 1 4 ) and C in (7.50a,d,e,f,h) is independent of α, as well as δ, h and ∆t. We are now in a position to prove the following convergence result concerning (P ∆t α,δ,h ). Theorem 7.2. Under all of the assumptions of Theorem 7.1, there exists a subsequence of {(u ∆t α,δ,h , σ ∆t α,δ,h , ̺ ∆t α,δ,h )} δ>0,h>0,∆t>0 , and functions

u α ∈ L ∞ (0, T ; H) ∩ L 2 (0, T ; V) ∩ W 1, 4 ϑ (0, T ; V ′ ) with u α (0) = u 0 , (7.51a) σ α ∈ L ∞ (0, T ; [L 2 (D)] 2×2 S ) ∩ L 2 (0, T ; [H 1 (D)] 2×2 S ) ∩ H 1 (0, T ; ([H 1 (D)] 2×2 S ) ′ ) (7.51b) with σ α non-negative definite a.e. in D T and σ α (0) = σ 0 , ̺ α ∈ L ∞ (0, T ; L 2 (D)) ∩ L 2 (0, T ; H 1 (D)) ∩ H 1 (0, T ; (H 1 (D)) ′ ) (7.51c) with ̺ α ≤ b a.e. in D T and ̺ α (0) = tr(σ 0 ), such that, as δ, h, ∆t → 0 + , u ∆t(,±) α,δ,h → u α weak* in L ∞ (0, T ; [L 2 (D)] 2 ), (7.52a) u ∆t(,±) α,δ,h → u α weakly in L 2 (0, T ; [H 1 (D)] 2 ), (7.52b) S ∂u ∆t α,δ,h ∂t → S ∂u α ∂t weakly in L 4 ϑ (0, T ; V), (7.52c) u ∆t(,±) α,δ,h → u α strongly in L 2 (0, T ; [L r (D)] 2 ), (7.52d) σ ∆t(,±) α,δ,h → σ α weak* in L ∞ (0, T ; [L 2 (D)] 2×2 ), (7.53a) σ ∆t(,±) α,δ,h → σ α weakly in L 2 (0, T ; [H 1 (D)] 2×2 ), (7.53b) E ∂σ ∆t α,δ,h ∂t → E ∂σ α ∂t weakly in L 2 (0, T ; [H 1 (D)] 2×2 S ), (7.53c) σ ∆t(,±) α,δ,h → σ α strongly in L 2 (0, T ; [L r (D)] 2×2 ), (7.53d) π h [β δ (σ ∆t(,±) α,δ,h )] → σ α strongly in L 2 (0, T ; [L 2 (D)] 2×2 ), (7.53e) Λ δ,m,p (σ ∆t(,±) α,δ,h ) → σ α δ mp strongly in L 2 (0, T ; [L 2 (D)] 2×2 ), m, p = 1, 2, (7.53f) and ̺ ∆t(,±) α,δ,h → ̺ α weak* in L ∞ (0, T ; L 2 (D)), (7.54a) ̺ ∆t(,±) α,δ,h → ̺ α weakly in L 2 (0, T ; H 1 (D)), (7.54b) E ∂̺ ∆t α,δ,h ∂t → E ∂̺ α ∂t weakly in L 2 (0, T ; H 1 (D)), (7.54c) ̺ ∆t(,±) α,δ,h → ̺ α strongly in L 2 (0, T ; L r (D)), (7.54d) π h β δ 1 - ̺ ∆t(,±) α,δ,h b → 1 - ̺ α b strongly in L 2 (0, T ; L 2 (D)), (7.54e) Λ δ,m,p 1 - ̺ ∆t(,±) α,δ,h b → 1 - ̺ α b δ mp strongly in L 2 (0, T ; L 2 (D)), m, p = 1, 2, (7.54f)
where ϑ ∈ (2, 4] and r ∈ [1, ∞).

Proof. The results (7.52a-c) follow immediately from the bounds (7.50a,f) on noting the notation (7.41a-c). The denseness of h>0 Q 1 h in L 2 (D) and (6.1d) yield that u α ∈ L 2 (0, T ; V). Hence, the result (7.51a) holds on noting (7.10) and (6.33), where u α : [0, T ] → H is weakly continuous.

The strong convergence result (7.52d) for u ∆t α,δ,h and r = 2 follows immediately from (7.50a) and the second bound in (7.50f) and (7.13) with X 0 = [H 1 (D)] 2 and X = X 1 = [L 2 (D)] 2 . Here we note that H 1 (D) is compactly embedded in L 2 (D). We note here also that one cannot appeal to (7.12) for this strong convergence result with

µ 0 = 2, µ 1 = 4/ϑ, Y 0 = [H 1 (D)] 2 , Y 1 = V ′ with norm S • H 1 (D) and Y = [L 2 (D)] 2 for the stated values of ϑ, since [L 2 (D)] 2 is not continuously embedded in V ′ as V is not dense in [L 2 (D)] 2 .
The result (7.52d) for u ∆t,± α,δ,h and r = 2 follows immediately from this result for u ∆t α,δ,h and the the bound on the last term on the left-hand side of (7.50a), which yields

u ∆t α,δ,h -u ∆t,± α,δ,h 2
L 2 (0,T ;L 2 (D)) ≤ C ∆t. (7.55)
Finally, we note from (7.11), for d = 2, that, for all η ∈ L 2 (0, T ; H 1 (D)), α,δ,h with r = 2 yields (7.52d) for u ∆t(,±) α,δ,h for the stated values of r. Similarly, the results (7.53a-c) follow immediately from (7.50b,g). The strong convergence result (7.53d) for σ ∆t α,δ,h follows immediately from (7.53b,h), (7.8) and (7.12) with µ

η L 2 (0,T ;L r (D)) ≤ C η 1-θ L 2 (0,T ;L 2 (D)) η θ L 2 (0,T ;H 1 (D)) (7.56) for all r ∈ [2, ∞) with θ = 1 -2 r ∈ (0,
0 = µ 1 = 2, Y 0 = [H 1 (D)] d×d , Y 1 = ([H 1 (D)] d×d ) ′ and Y = [L r (D)
] d×d for the stated values of r. Here we note that H 1 (D) is compactly embedded in L r (D) for the stated values of r, and L r (D) is continuously embedded in (H 1 (D)) ′ . Similarly to (7.55) and (7.56), the last bound in (7.50b) then yields that (7.53d) holds for σ ∆t(,±) α,δ,h . The results (7.54a-d) follow analogously from noting (7.50c,g). Hence, on noting (7.8) and (6.33), the results (7.51b,c) hold, where σ α : [0, T ] → [L 2 (D)] d×d S and ̺ α : [0, T ] → L 2 (D) are weakly continuous, apart from the claims on the non-negative definiteness of σ α and the bound on ̺ α . It remains to prove these, (7.53e,f) and (7.54e,f). It follows from (1.11), (6.52b) and (7.50b,e) that

[σ α ] -L 2 (0,T ;L 1 (D)) (7.57) ≤ [σ α ] --[σ ∆t(,±) α,δ,h ] - L 2 (0,T ;L 1 (D)) + [σ ∆t(,±) α,δ,h ] --π h [σ ∆t(,±) α,δ,h ] - L 2 (0,T ;L 1 (D)) + π h [σ ∆t(,±) α,δ,h ] - L 2 (0,T ;L 1 (D)) ≤ σ α -σ ∆t(,±) α,δ,h L 2 (0,T ;L 1 (D)) + C [h + δ] .
The desired non-negative definiteness result on σ α in (7.51b) then follows from (7.53d). The desired bound on ̺ α in (7.51c) follows similarly from (1.11), (6.52b), (7.50c,e) and (7.54d). The results (7.53e,f) follow immediately from (6.52a), (7.50b), (1.11), (7.53d), (3.3) and the nonnegative definiteness result on σ α in (7.51b). The results (7.54e,f) follow similarly from the scalar version of (6.52a), (7.50c), (1.11), (7.54d), (3.3) and the bound on ̺ α in (7.51c).

Lemma 7.3. Under all of the assumptions of Theorem 7.1, the subsequence of {(σ ∆t α,δ,h , ̺ ∆t α,δ,h )} δ>0,h>0,∆t>0 of Theorem 7.2 and the limiting functions σ α and ̺ α , satisfying (7.51b,c), are such that, as δ, h, ∆t → 0 + ,

π h β b δ (̺ ∆t(,±) α,δ,h ) → ̺ α = tr(σ α ) strongly in L 2 (0, T ; L 2 (D)). (7.58)
In addition, we have with h = o(δ), as δ → 0 + , that σ α is positive definite and tr(σ α ) < b a.e. in D T . (7.59) Proof. Choosing φ = η I in (7.43b) and subtracting from (7.43c) yields that

DT π h ∂ ∂t ̺ ∆t α,δ,h -tr(σ ∆t α,δ,h ) η dx dt + α DT ∇ ̺ ∆t,+ α,δ,h -tr(σ ∆t,+ α,δ,h ) • ∇η dx dt (7.60) + DT d m=1 d p=1 [u ∆t,- α,δ,h ] m b Λ δ,m,p 1 - ̺ ∆t,+ α,δ,h b + tr Λ δ,m,p (σ ∆t,+ α,δ,h ) ∂η ∂x p dx dt = 0 ∀η ∈ L 2 (0, T ; Q 1 h ).
It follows from (7.52d), (7.53b-d,f), (7.54b-e), (7.7), (7.51a-c) and (6.9a,b) that we may pass to the limit δ, h, ∆t → 0 + in (7.60) with η = π h χ to obtain

T 0 ∂ ∂t (̺ α -tr(σ α )) , χ H 1 (D) dt + α DT ∇ (̺ α -tr(σ α )) • ∇χ dx dt (7.61) - DT (̺ α -tr(σ α )) u α • ∇χ dx dt = 0 ∀χ ∈ C ∞ 0 (0, T ; C ∞ (D)),
where [̺ α -tr(σ α )](0) = 0. For example, in order to pass to the limit on the first term in (7.60), we note that

DT π h ∂ ∂t ̺ ∆t α,δ,h -tr(σ ∆t α,δ,h ) π h χ dx dt (7.62) = DT ∂ ∂t ̺ ∆t α,δ,h -tr(σ ∆t α,δ,h ) π h χ + (I -π h ) ̺ ∆t α,δ,h -tr(σ ∆t α,δ,h ) π h ∂χ ∂t dx dt.
Hence the desired first term in (7.61) follows from noting (7.53b,c), (7.54b,c), (7.7) and (6.9a,b). As C ∞ 0 (0, T ; C ∞ (D)) is dense in L 2 (0, T ; H 1 (D)), we have, on noting (7.51a-c), that (7.61) holds for all χ ∈ L 2 (0, T ; H 1 (D)). It then follows from (7.51a-c) that we can choose χ = ̺ α -tr(σ α ) in (7.61) to yield that ̺ α = tr(σ α ) as [̺ α -tr(σ α )](0) = 0. Recalling (7.51b,c), we have that ̺ α ∈ [0, b] a.e. in D T . The desired convergence result (7.58) then follows from noting this, (5.3), (3.3), (6.52b), (7.50c) and (7.54d).

We now improve on (

π h tr(β δ (σ ∆t,+ α,δ,h )) dx dt + D b T tr σ α -π h β δ (σ ∆t,+ α,δ,h ) dx dt =: T 1 + T 2 .
We deduce from (3.6) and (3.2) that

T 1 = D b T π h tr A δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) β δ (σ ∆t,+ α,δ,h ) + 2 β δ 1 - ̺ ∆t,+ α,δ,h b dx dt. (7.64)
It follows from (1.3b) and (1.4b) that tr A δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) β δ (σ ∆t,+ α,δ,h ) ≤ tr (A δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h )) 2 β δ (σ ∆t,+ α,δ,h ) tr β δ (σ ∆t,+ α,δ,h ) L 2 (0,T ;H 1 (D)) .

It follows from (7.69), (7.51c), (6.26), (3.3), (6.8), (7.50c) and (7.54e) that T 1 = 0. In addition, it follows immediately from (7.53e) that T 2 = 0. Hence, we conclude from (7.63) that |D b T | = 0, and so ρ α = tr(σ α ) < b a.e. in D T ; that is, the second desired result in (7.59).

We now establish the other result in (7.59) that σ α is symmetric positive definite a.e. in D T , which improves on (7.51b). This result requires the further assumption that h = o(δ), as δ → 0 + . Assume that σ α is not symmetric positive definite a. A simple variation of (6.7a), (1.3b), (1.4a,b) and (7.50d) yield that T 3 ≤ DT π h A δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) [β δ (σ ∆t,+ α,δ,h )] , where we have noted from (1.3b), (1.4b), (1.8), (6.51), (3.6) ≤ C.

We will now show, on possibly extracting a further subsequence of {(σ ∆t α,δ,h , ̺ ∆t α,δ,h )} δ>0,h>0,∆t>0 , that

π h G ′ δ 1 - ̺ ∆t,+ α,δ,h b β δ (σ ∆t,+ α,δ,h ) → 1 - ρ α b -1
σ α weakly in L ) to some limit for a subsequence. We just need to show it is the limit stated in (7.75). We have from (7.53e), (7.54d Similarly to (7.72), we have that (7.73), (7.75) and our definitions of χ ∆t,+ δ,h and v yield that T 4 = 0. Hence it follows from (7.70) with T 3 = T 4 = 0 that |D 0 T | = 0, and so σ α is positive definite a.e. in D T ; that is, the first desired result in (7.59). Lemma 7.4. Under all of the assumptions of Lemma 7.3, a further subsequence of the subsequence of {(σ ∆t α,δ,h , ̺ ∆t α,δ,h )} δ>0,h>0,∆t>0 of Lemma 7.3 and the limiting function σ α , satisfying (7.51b) and (7.59), are such that, as δ, h, ∆t → 0 + , with h = o(δ), π h κ δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) β δ (σ ∆t,+ α,δ,h ) → σ α strongly in L 2 (0, T ; [L 2 (D)] 2×2 ), (7.77a) π h A δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) β δ (σ ∆t,+ α,δ,h ) → A(σ α ) σ α weakly in L ), (7.77b) π h κ δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) A δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) β δ (σ ∆t,+ α,δ,h ) → A(σ α ) σ α (7.77c) weakly in L 2 (0, T ; [L 2 (D)] 2×2 ).

Proof. It follows from (6.7a), (1.3b), (1.4b), (5.2), (1.8), (3.3), (6.8) and (7.50b) that π h κ δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) β δ (σ ∆t,+ α,δ,h ) -π h β δ (σ ∆t,+ α,δ,h ) × π h A δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) β δ (σ ∆t,+ α,δ,h ) .

It follows immediately from (7.79), (7.53e), (7.58) and (7.50h) that the weak limits in (7.77b,c) are the same. Hence, the desired result (7.77c). Proof. The function spaces and the initial conditions for (u α , σ α ) follow immediately from (7.51a,b), (7.59) and (7.77c). It remains to prove that (u α , σ α ) satisfy (7.80a,b). It follows from (6.3), (7.52b-d), (7.77c), (4.11b), (7.9) and (4.14) that we may pass to the limit, δ, h, ∆t → 0 + , with h = o(δ), in (7.43a) to obtain that (u α , σ α ) satisfy (7.80a).

It follows from (7.53b-f), (7.52b,d), (7.77a,b), (7.7), (6.9a,b), (1.5a) and as u α ∈ L 2 (0, T ; V) that we may pass to the limit δ, h, ∆t → 0 + , with h = o(δ), in (7.43b) with χ = π h φ to obtain (7.80b) for any φ ∈ C ∞ 0 (0, T ; [C ∞ (D)] 2×2 S ). For example, similarly to (7.62), in order to pass to the limit on the first term in (7.43b) ). Of course passing to the limit δ, h, ∆t → 0 + , with h = o(δ), in (7.43c), using in addition (7.54b-d,f), yields the weak formulation for tr(σ α ) consistent with (7.80b).

Remark 7.2. Choosing 1 2 (φ + φ T ) as a test function in (7.80b) for any φ ∈ L 2 (0, T ; [H 1 (D)] 2×2 ) yields, on noting the symmetry of σ α , (7.80b) with the term 2 (∇u α )σ α replaced by (∇u α )σ α + σ α (∇u α ) T , which is consistent with (5.1c).

Finally, it follows from (7.50a,f,h), (7.52a-c) and (7.77c) that

1 . Introduction 1 . 1 . 1 I

 1111 The FENE-P model. We consider the standard FENE-P model for a dilute polymeric fluid. The fluid, confined to an open bounded domain D ⊂ R d (d = 2 or 3) with a Lipschitz boundary ∂D, is governed by the following non-dimensionalized system for a given b ∈ R >0 :(P) Find u : (t, x) ∈ [0, T ) × D → u(t, x) ∈ R d , p : (t, x) ∈ D T := (0, T ) × D → p(t, x) ∈ R and σ : (t, x) ∈ [0, T ) × D → σ(t, x) ∈ R d×d S,>0,b such that Re ∂u ∂t + (u • ∇)u = -∇p + (1 -ε)∆u + ε Wi div (A(σ) σ) + f on D T , (1.1a) div u = 0 on D T , (1.1b) ∂σ ∂t + (u • ∇)σ = (∇u)σ + σ(∇u) T -A(σ) σ Wi on D T , (1.1c) u(0, x) = u 0 (x) ∀x ∈ D, (1.1d) σ(0, x) = σ 0 (x) φ -1 ∀φ ∈ R d×d S,>0,b := ψ ∈ R d×d S,>0 : tr(ψ) < b . (1.2)

  • [H 1 (D)] d and | • | [H 1 (D)] d , respectively. We denote the duality pairing between H -1 (D) and H 1 0 (D) as •, • H 1 0 (D) , and we similarly write •, • H 1 0 (D) for the duality pairing between e.g. [H -1 (D)] d and [H 1 0 (D)] d . For notational convenience, we introduce also convex sets such as [H 1 (D)] d×d S,>0 := {φ ∈ [H 1 (D)] d×d S : φ ∈ R d×d S,>0 a.e. in D}, and [H 1 (D)] d×d S,>0,b := {φ ∈ [H 1 (D)] d×d S : φ ∈ R d×d S,>0,b a.e. in D}.

Corollary 2 . 1 .

 21 Under the assumptions of Proposition 2.1 it follows that

  sup

Corollary 3 . 1 .

 31 Under the assumptions of Proposition 3.1 it follows that sup t∈(0,T )

(4. 21 )

 21 Combining (4.18)-(4.21) yields the first desired inequality in(4.16

1 2 .

 2 Then (3.15), (3.16), (3.5), (4.30) and (4.29) yield that

  .50) but this contradicts (4.48) on recalling (3.2) and (3.3). Therefore, tr(σ n h )| K k < b on all simplices K k , and so it follows from (4.45), (3.2) and (3.3) that

δ b 2 aRemark 5 . 1 .

 251 .e. in D T .Combining (5.9)-(5.13a,b) yields the desired result(5.8). We note if we multiply the advection term in (5.5c) by -G ′ δ (1 -tr(σ α,δ ) b

  ([H 1 (D)] d×d S ) ′ be the topological dual of [H 1 (D)] d×d S with [L 2 (D)] d×d S being the pivot space. Let E : ([H 1 (D)] d×d S ) ′ → [H 1 (D)] d×d S be such that Eχ is the unique solution of the Helmholtz problem D [∇(Eχ) :: ∇φ + (Eχ) : φ] dx = χ, φ H 1 (D) ∀φ ∈ [H 1 (D)] d×d S , (7.7) where •, • H 1 (D) denotes the duality pairing between ([H 1 (D)] d×d S ) ′ and [H 1 (D)] d×d

  39) by ∆t n , summing from n = 1, ..., N T and noting (1.3b), (1.4b), (6.48), (1.8), (3.3), (6.51),(7.11) with d = 2, (6.8) and (7.15a) yields that NT n=1 ∆t n π h A δ (σ n α,δ,h , ̺ n α,δ,h ) β δ (σ n α,δ,h )

4 ϑL 2 (

 42 D) dt ≤ C,(7.49)

  7.51c) by establishing that ̺ α = tr(σ α ) < b a.e. in D T . Assuming that tr(σ α ) = b a.e. in D b T ⊂ D T , we have b |D b T | =

 8 . 4 L 4 ( 2 L 2 ( 2 H 1 (T π h β δ 1 -dx dt ≤ C π h β δ 1 -

 844222111 64) and (7.65) yields, on noting a scalar version of (6.7a) over D b T , (7.50d),(1tr β δ (σ ∆t,+ α,δ,h ) + 1 β δ 1 -̺ ∆t,+ α,δ,hIt follows from (6.51),(7.11), for d = 2, and (7.50b) thatDT π h σ ∆t,+ α,δ,h 4 dx dt ≤ C σ ∆t,+ α,δ,h D) σ ∆t,+ α,δ,h D) dt ≤ C. (7.67)Similarly to (6.51), it follows from (6.4) withK k replaced by (K k × (t n-1 , t n ))∩D b T , k = 1, . . . , N K and n = 1, . . . , N T , that D b where z := π h β δ 1 -̺ ∆t,+ α,δ,h b -1 -̺α b .Here we have noted that ̺ α = tr(σ α ) = b a.e. in D b T . Combining (7.66)-(7.68) yields, on noting(7.11), for d = 2, thatT 1 ≤ C z

2 3 H 1 (D) dt 3 8 ≤ C z 3 4 L 3 4 L 2 ( 2 L 2 4 L

 23184342224 (0,T ;L 2 (D)) z 1 (0,T ;L 2 (D)) z 1 ∞ (0,T ;L 2 (D)) z 1 4

DT π h G ′ δ 1 -

 1 e. in D 0 T ⊂ D T . Let v ∈ L ∞ (0, T ; [L ∞ (D)] 2 ) be such that σ α v = 0 with v = 1 a.e. in D 0 T and v = 0 a.e. in D T \ D 0 T . It then follows from (3.2) and (3.6) that|D 0 T | = DT v dx dt = DT π h G ′ δ (σ ∆t,+ α,δ,h ) β δ (σ ∆t,+ α,δ,h ) v dx dt (7.70) ≤ DT π h A δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) β δ (σ ∆t,+ α,δ,h ) v dx dt + ̺ ∆t,+ α,δ,h b β δ (σ ∆t,+ α,δ,h ) v dx dt =: T 3 + T 4 .

π

  h β δ (σ ∆t,+ α,δ,h ) :: (v v T ) δ (σ ∆t,+ α,δ,h ) :: (v v T ) dx dt

1 2 .π

 12 Then(7.53e) and the definition of v yield that, as δ, h, ∆t → 0 + ,DT π h β δ (σ ∆t,+ α,δ,h ) :: (v v T ) dx dt → DT σ α :: (v v T) dx dt = 0, (7.72) so we have that T 3 = 0. Similarly to (7.71), on setting χ ∆t,+ h β δ (σ ∆t,+ α,δ,h ) , we have from (6.7a) thatT 4 ≤ DT π h [χ ∆t,+ δ,h ]

8 5 8 5

 88 (0, T ; [L (D)] 2×2 ). (7.75) as δ, h, ∆t → 0 + with h = o(δ). It follows immediately from (7.74) and our definition of χ ∆t,+ δ,h that π h χ ∆t,+

  ) and (3.1) that χ ∆t,+ δ,h → 1 -ρα b -1 σ α a.e. on D T , for a subsequence, as we have already established that ρ α = tr(σ α ) < b a.e. on D T . So it remains to establish that (I -π h ) χ ∆t,+ δ,h converges to zero a.e. on D T . As G ′δ ∈ C 0,1 (R) is monotonic, it follows from (6.9b), (6.53), (6.4), (3.2) and (7.50b,c) that(I -π h ) χ ∆t,+ δ,h L 1 (DT ) ≤ (I -π h ) π h G ′ δ 1 -̺ ∆t,+ α,δ,h b π h β δ (σ ∆t,+ α,δ,h ) L 1 (DT ) (7.76) + (I -π h ) G ′ δ 1 -̺ ∆t,+ α,δ,h b L 2 (DT ) π h β δ (σ ∆t,+ α,δ,h ) L 2 (DT ) ≤ C h ∇π h G ′ δ 1 -̺ ∆t,+ α,δ,h b L 2 (DT ) π h β δ (σ ∆t,+ α,δ,h ) L 2 (DT ) ≤ C δ -1 h.Hence, we have for a subsequence that (I -π h ) χ ∆t,+ δ,h converges to zero a.e. on D T as δ, h, ∆t → 0 + with h = o(δ). Therefore we have established(7.75).

≤ 2 L 1 2 L 2

 2122 C π h β b δ (̺ ∆t,+ α,δ,h ) -tr(β δ (σ ∆t,+ α,δ,h ) 1 (DT ) ≤ C π h β b δ (̺ ∆t,+ α,δ,h ) -tr(β δ (σ ∆t,+ α,δ,h ) 1 (DT )

Theorem 7 . 3 . 1 -L 4 4DT 2 (

 73142 Under all of the assumptions of Lemma 7.3, the limiting functions (u α , σ α ) satisfying (7.51a,b) and (7.59) solve the following problem:(P α ) Find u α ∈ L ∞ (0, T ; H) ∩ L 2 (0, T ; V) ∩ W 1, 4 ϑ (0, T ; V ′ ) and σ α ∈ L ∞ (0, T ; [L 2 (D)] 2×2 S,>0,b ) ∩ L 2 (0, T ; [H 1 (D)] 2×2 S ) ∩ H 1 (0, T ; ([H 1 (D)] 2×2 S ) ′ ), with A(σ α ) σ α ∈ L 2 (0, T ; ([L 2 (D)] 2×2 S )), such that u α (0) = u 0 , σ α (0) = σ 0 and Re ε) ∇u α : ∇v + Re [(u α • ∇)u α ] • v] dx dt (α ) σ α : ∇v dx dt ∀v ∈ α • ∇)σ α : φ + α ∇σ α :: ∇φ] dx dt (7.80b) = ∇u α ) σ α -1 Wi A(σ α ) σ α : φ dx dt ∀φ ∈ L 2 (0, T ; [H 1 (D)] 2×2 S ),where ϑ ∈ (2, 4).

+ A(σ α ) σ α 2 L 2 (

 22 D) dt ≤ C,(7.82) where ϑ ∈ (2, 4) and C is independent of the stress diffusion coefficient α. Of course, in addition, it follows from (7.59) and (1.8) that σ α L ∞ (DT ) < b.

  3.35a,b) there.

	Theorem 4.2. For all regular partitionings T h of D into simplices {K k } NK k=1 and all partitionings {∆t n } NT n=1 of [0, T ], there exists a subsequence {{(u n δ,h , σ n δ,h )} NT n=1 } δ>0 , where {(u n δ,h , σ n δ,h )} NT n=1 ∈ [V 0 h × S 0

h ] NT solves (P ∆t δ,h ), (4.13a,b), and

  .24) Combining (7.22)-(7.24) and (7.19), yields, on applying a Young's inequality, that for all ζ > 0

  and (7.50i) that≤ C + C π h A δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) β δ (σ ∆t,+ α,δ,h )

	(7.74)	DT	π h [χ ∆t,+ δ,h ]	1 2	2	dx dt ≤ C π h χ ∆t,+ δ,h	L 1 (DT )
						≤ C π h χ ∆t,+ δ,h	L 1 (DT )	≤ C π h χ ∆t,+ δ,h	L	8 5 (DT )
								L	8 5 (DT )

  (κ δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) -1) 2 tr(β δ (σ ∆t,+ α,δ,h )) π h tr(β δ (σ ∆t,+ α,δ,h )) dx dtThe desired result (7.77a) then follows immediately from (7.78), (7.53e) and (7.58).The desired result (7.77b) follows immediately from (3.6), (3.2), (7.75) and (1.2) as ρ α = tr(σ α ).Similarly to (7.78), it follows from (6.7a), (1.3b), (1.4b), (5.2), (6.8), and (7.50d) thatπ h κ δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) A δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) β δ (σ ∆t,+ α,δ,h ) -π h A δ (σ ∆t,+α,δ,h , ̺ ∆t,+ α,δ,h ) β δ (σ ∆t,+ α,δ,h )

	(7.78)			2 L 2 (DT )
	≤ π h ≤ C DT DT	π h β b δ (̺ ∆t,+ α,δ,h ) -tr(β δ (σ ∆t,+ α,δ,h )	π h σ ∆t,+ α,δ,h	+ δ dx dt
	≤ C π h β b δ (̺ ∆t,+ α,δ,h ) -tr(β δ (σ ∆t,+ α,δ,h )	L 2 (DT )	.
	(7.79)
					L 1 (DT )
	≤ C

DT π h β b δ (̺ ∆t,+ α,δ,h ) -tr(β δ (σ ∆t,+ α,δ,h )

1 2

  , we note that The desired result (7.80b) then follows from noting thatC ∞ 0 (0, T ; [C ∞ (D)] 2×2 S ) is dense in L 2 (0, T ; [H 1 (D)] 2×2 S

	(7.81)
	DT	π h	∂σ ∆t α,δ,h ∂t	+	A δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) β δ (σ ∆t,+ α,δ,h ) Wi	: π h φ dx dt
		=	DT	 	∂σ ∆t α,δ,h ∂t	+	π h A δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) β δ (σ ∆t,+ α,δ,h ) Wi	  : π h φ dx dt
			+	DT	(I -π h )	  σ ∆t α,δ,h : π h	∂φ ∂t	-	π h A δ (σ ∆t,+ α,δ,h , ̺ ∆t,+ α,δ,h ) β δ (σ ∆t,+ α,δ,h ) Wi	 : π h φ  dx dt.