# New Schemes, Positive Definiteness and Energy Preserving Approaches for Oldroyd-B Model Under Slip Boundary Condition

J K Djoko, J Koko, T Sayah

► To cite this version:

J K Djoko, J Koko, T Sayah. New Schemes, Positive Definiteness and Energy Preserving Approaches for Oldroyd-B Model Under Slip Boundary Condition. International Journal of Applied Mathematics, Computational Science and Systems Engineering, 2025, 7, pp.1 - 18. 10.37394/232026.2025.7.1. hal-05055430

HAL Id: hal-05055430

https://hal.science/hal-05055430v1

Submitted on 5 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# New Schemes, Positive Definiteness and Energy Preserving Approaches for Oldroyd-B Model Under Slip Boundary Condition

J.K. DJOKO¹, J. KOKO², T. SAYAH³

¹School of Mathematical and Statistical Sciences,
North West University
SOUTH AFRICA

²LIMOS, Université Clermont-Auvergne - CNRS UMR 6158,
B.P 10448, 63000 Clermont-Ferrand,
FRANCE

³Laboratoire de Mathematiques et Applications,
Unite de recherche "Mathematiques et Modelisation",
Faculté des Sciences, Université Saint-Joseph, B.P 11-514 Riad El Solh,
Beirut 1107 2050,
LEBANON

*Abstract:* New numerical schemes for simulating the Oldroyd-B model are formulated and analyzed. The first approach is based on the semi-implicit Euler scheme, while the second method applies the Marchuk-Yanenko splitting scheme. The goal is to develop schemes that preserve the symmetry and positive definiteness of the conformation tensor, as well as the energy estimate, which includes kinetic energy and the elastic energy of the extra stress, as previously obtained by A. Lozinski and R. Owens in J. Non-Newtonian Fluid Mech., **112**, 161–176, 2003. We finally present numerical results which are consistent with theoretical predictions.

*Key-Words:* Navier slip boundary condition, Oldroyd model, energy estimate, semi-implicit scheme, Marchuk-Yanenko.

AMS subject classification: 65M12, 65N30, 76M10, 35B45, 76A10, 35B35.

Received: March 2, 2024. Revised: August 12, 2024. Accepted: December 3, 2024. Published: March 12, 2025.

## 1 Introduction

Let $\Omega \subset \mathbb{R}^d$ (with $d = 2,3$) be an open, bounded domain with Lipschitz-continuous boundary $\Gamma = \partial\Omega$ and let $(0, T)$ be a time interval. We consider a viscoelastic fluid of Oldroyd type and assume that the body force acting on the fluid is given by $b$ per unit mass. Without loss of generality, we assume that the fluid's mass density is equal to unity. The goal is to determine the velocity $u(x, t)$, the pressure $p(x, t)$, and the elastic stress $E(x, t)$, which satisfy the following set of equations:

$$ 
\begin{aligned}
& \mathrm{Re}(u_t + (u \cdot \nabla)u) - 2\beta \operatorname{div} Du + \nabla p \\
& \qquad = \operatorname{div} E + b \quad \text{in } \Omega \times (0, T), \\
& \operatorname{div} u = 0 \quad \text{in } \Omega \times (0, T),
\end{aligned}
$$

$$ 
\begin{aligned}
& \frac{1}{We}E + E_t + (u \cdot \nabla)E - \nabla uE - E(\nabla u)^T \\
& \qquad = 2\frac{1-\beta}{We}Du \quad \text{in } \Omega \times (0, T),
\end{aligned}
$$

$$ E = E^T, \qquad (1.1) $$

where $\beta \in (0, 1)$ is the dimensionless solvent parameter (regarded as the total viscosity ratio),

and $\mathrm{Re}$ and $\mathrm{We}$ are positive physical parameters representing the Reynolds number and the Weissenberg number. The symmetric part the velocity gradient is $2Du = \nabla u + (\nabla u)^T$. The Oldroyd-B constitutive model, as expressed in the third equation of (1.1), represents the most straightforward nonlinear extension of Maxwell's concept, which involves defining stress through a system of ordinary differential equations based on the velocity gradient and its time derivative. This fundamental model for complex fluids highlights the significant challenge posed by high Weissenberg number regimes. The formulation of the third equation in (1.1) aims to establish a relationship between the stress tensor $E$ and the strain rate tensor $Du$. The Weissenberg number serves as a key parameter that differentiates viscoelastic fluids from Newtonian fluids. To solve (1.1), appropriate boundary and initial conditions must be specified. In this context, we assume that the boundary $\Gamma$ of $\Omega$ is impermeable. Hence

$$ u \cdot n = 0 \quad \text{on } \Gamma \times (0, T), \qquad (1.2) $$

where $n$ is the outward unit normal vector to $\Gamma$. The Cauchy stress tensor within the fluid is $T$ given by
$T = -pI + 2\beta D u + E$, with $I$ the d-dimensional identity matrix. Next, we describe the relation between the stress and the velocity on $\Gamma$. Taking the scalar product of $u$ with the first equation in (1.1), we obtain

$$ \frac{Re}{2} \frac{d}{dt} \int_{\Omega} |u|^2 dx + \int_{\Omega} T : Du dx \\ + \int_{\Gamma} (-Tn)_{\tau} \cdot u_{\tau} d\sigma = \int_{\Omega} b \cdot u dx \quad (1.3) $$

where $d\sigma$ is the surface measure associated to $\Gamma$. We recall that for any vector $w$ defined on $\Gamma$, we set $w_{\tau} = w - (w \cdot n)n$. Thus $(Tn)_{\tau}$ denotes the projection of the normal stress onto the corresponding tangent plane. In (1.3), $\frac{Re}{2} \frac{d}{dt} \int_{\Omega} |u|^2 dx$ represents the variation of kinetic energy, $\int_{\Omega} T : Du dx$ is the dissipation mechanisms in bulk, while $\int_{\Gamma} (-Tn)_{\tau} \cdot u_{\tau} d\sigma$ stands for dissipative processes on the $\Gamma$. The most general relation between $u_{\tau}$ and $(Tn)_{\tau}$ is to consider the implicit constitutive relation (see [33])

$$ \psi(u_{\tau}, (Tn)_{\tau}) = 0 $$

where $\psi$ is a function. If $(Tn)_{\tau} = 0$, this corresponds to a perfect slip boundary condition, whereas if $u_{\tau} = 0$, there is no slip. In this work, we assume that the liquid-solid interaction is governed by Navier's slip condition

$$ (tn)_{\tau} + \alpha(u_{\tau} - w_{\tau}) = 0 \quad \text{on } \Gamma \times (0, T), \quad (1.4) $$

where $\alpha$ is a positive function and $w$ is the velocity of the solid surface satisfying $w \cdot n|_{\Gamma} = 0$. It is worth observing that $w_{\tau} \neq 0$ indicate the fact that the flow region $\Omega$ is fixed but its boundary $\Gamma$ is not; in fact, it may undergo tangential motion. Finally, we assume that

$$ u(x, 0) = u_0, \quad \text{and} \quad E(x, 0) = E_0, \quad \text{in } \bar{\Omega} \quad (1.5) $$

where $u_0$ and $E_0$ are given functions whose regularity will be specified later, with $\text{div } u_0 = 0$.

Many research works have explored the numerical aspects (such as simulations and the convergence of numerical schemes) as well as the mathematical properties (including existence, non-existence, and regularity) of viscoelastic fluids modeled by the Oldroyd-B equations in two or three dimensions, under Dirichlet or periodic boundary conditions. (see [1, 2, 3, 4, 5, 6, 8, 9, 14, 15, 17, 18, 19, 20, 21, 22, 23, 26, 27, 28, 29, 30, 31, 32, 34, 38, 39]). These studies focus on time-dependent or steady

flows of incompressible fluids in bounded domains, under various restrictions on the material constants and other data. In [35], the study focuses on the existence of regular solutions for the steady flow of an Oldroyd-type viscoelastic fluid, assuming slip conditions on the entire boundary.

In [19], the authors investigate an Oldroyd-B model with Dirichlet boundary conditions, proposing a scheme that preserves the quantity $Re \|u\|^2 + \int_{\Omega} \text{tr} E dx$. Numerical examples are provided to showcase the advantages of their method over classical approximations. Additionally, in [13], we develop a numerical scheme that preserves symmetry and positive definiteness, and provide an a priori estimate for the Oldroyd-B model with Tresca boundary conditions. This contribution goes beyond a straightforward extension of [19], as the non-standard boundary condition leads to a mixed variational inequality of the second kind. In this work,

In this contribution, we considered the Oldroyd-B model with slip boundary condition of Navier's type. We begin by establishing a basic energy estimate. Next, we formulate two new discretizations: the first scheme is based on semi-implicit Euler's method while the second method is based on Marchuk-Yanenko's decomposition algorithm. The semi-implicit scheme we formulate linearizes the fully implicit Euler scheme presented in [19] and leads directly to a numerical algorithm based on a space and time discretization. It is important to highlight that the presented semi-implicit scheme is consistent, while the decay of free energy serves as a form of stability result. The second scheme we propose in this work is a direct application of the Marchuk-Yanenko splitting approach (see [23, 32, 24]). This new scheme is designed to ensure free energy dissipation at the discrete level (see Proposition 4.1). To our knowledge, this approach has not been previously applied—at least in this context—to assess the numerical stability of the Oldroyd-B model. The key idea is to formulate sub-problems that not only lead to the decay of free energy but also simplify the computation of the unknowns. The implementations are carried out using the FreeFem code [25].

The remainder of this work is organized as follows:

* In Section 2, we derive the a priori estimate.

* Section 3 focuses on the space and time discretization of the problem. We prove the discrete counterpart of the energy estimate.

* Section 4 addresses the second time discretization of the problem, based on Marchuk-Yanenko's algorithm and proves the discrete counterpart of the
energy estimate.

• Section 5 is dedicated to the numerical validation of the theoretical results.

## 2 A priori Analysis

To the best of our knowledge, the mathematical analysis of (1.1)-(1.5) has not yet been thoroughly studied. Therefore, in the absence of a solid theoretical framework to develop efficient and reliable numerical solutions, it is hard to make strong statements about the numerical results. However, we will present a basic energy estimate in the spirit of [19, 8] for the continuous equations. Our aim is to derive certain a priori estimates for the solution. The constitutive equation for the extra stress tensor E can be expressed in terms of the conformation tensor

$$A = \frac{We}{1 - \beta} E + I; \text{ with } A \text{ given by the relation}$$

$$A_t + (\mathbf{u} \cdot \nabla)A - \nabla \mathbf{u}A - A(\nabla \mathbf{u})^T = -\frac{1}{We}A + \frac{1}{We}I. \quad (2.1)$$

(2.1) is the starting point for most derivations of logarithmic reformulations, and other analyses. It has been established in [8] that

**Lemma 2.1.** If $A(t=0) = \frac{We}{1-\beta}E(t=0) + I$ is positive definite and symmetric, then this property is propagated forward in time by (2.1).

The theoretical analysis of time dependent problems usually rely on the following Gronwall's lemma (see for instance [36] and [12])

**Lemma 2.2.** Let $T > 0$, and $\phi$ be a non-negative function in $L^1(0, T)$. Let $c$ be a positive constant and $\phi \in C([0, T])$ a function satisfying

$$\forall t \in [0, T] \quad 0 \le \phi(t) \le c + \int_0^t \alpha(s)\phi(s)ds.$$

Then $\phi$ verifies the bound

$$\forall t \in [0, T] \quad \phi(t) \le c \exp \left( \int_0^t \alpha(s)ds \right).$$

We have the following a priori estimate

**Proposition 2.1.** Let $b \in L^\infty(0, \mathbb{R}; L^2(\Omega)^d)$, $\alpha \in L^\infty(\Gamma)$ and $w_\tau \in L^2(\Gamma)$. Let $(\mathbf{u}, p, A)$ be a regular solution of $(1.1)_1$, $(1.1)_2$, $(1.2)$, $(1.4)$ and (2.1). Assume that there exists $\delta$ such that $\beta < \delta < \alpha$.

Then there is a constant $c$ depending on $\Omega$ such that

$$\begin{aligned} & \mathrm{Re}\|\mathbf{u}\|^2 + \frac{1-\beta}{We} \int_{\Omega} \mathrm{tr} A \\ & \le \exp(-\gamma t) \left[ \mathrm{Re}\|\mathbf{u}_0\|^2 + \frac{1-\beta}{We} \int_{\Omega} \mathrm{tr} A_0 \right] \\ & \quad + \frac{1-\beta}{We^2} d|\Omega| \frac{1}{\gamma} (1 - \exp(-\gamma t)) \\ & \quad + \frac{1}{\gamma} \|\alpha\|_{L^\infty(\Gamma)} \|w_\tau\|_{L^\infty(0,t;L^2(\Gamma))}^2 (1 - \exp(-\gamma t)) \\ & \quad + \frac{c}{\gamma\beta} (1 - \exp(-\gamma t)) \|b\|_{L^\infty(0,t;L^2(\Omega))}^2, \end{aligned}$$

$$\text{with } \gamma = \min \left( \frac{c}{\mathrm{Re}}, \frac{1}{We} \right).$$

**proof:** We introduce the following double contraction between rank-two tensors $A, B \in \mathbb{R}^{d \times d}$:

$$A:B = \mathrm{tr}(AB^T) = \mathrm{tr}(A^T B) = \sum_{1 \le i,j \le d} A_{ij} B_{ij}.$$

We define the following spaces:

$$\mathcal{V}_{\tau} = \{ \mathbf{v} \in \mathcal{H}^1(\Omega) : \mathbf{v} \cdot \mathbf{n}|_{\Gamma} = 0 \}$$

and

$$M = {}_0^2(\Omega) = \{q \in L^2(\Omega), (q, 1) = 0\}.$$

Let $(\mathbf{v}, q) \in \mathcal{V}_{\tau} \times L_0^2(\Omega)$ and $\psi \in \{\psi \in H^1(\Omega)^{d \times d}, \psi = \psi^T\}$. Then one verifies easily that $(\mathbf{u}, p, A)$ solves

$$\begin{aligned} & \mathrm{Re}(\mathbf{u}_t, \mathbf{v}) + \mathrm{Re}((\mathbf{u} \cdot \nabla)\mathbf{u}, \mathbf{v}) + 2\beta(D\mathbf{u}, D\mathbf{v}) \\ & \quad - (p, \mathrm{div}\,\mathbf{v}) = -\frac{1-\beta}{We}(A, D\mathbf{v}) \\ & \quad - \int_{\Gamma} \alpha(\mathbf{u}_{\tau} - \mathbf{w}_{\tau}) \cdot \mathbf{v}_{\tau} + \int_{\Omega} \mathbf{b} \cdot \mathbf{v} \\ & (q, \mathrm{div}\,\mathbf{u}) = 0 \\ & (A_t, \psi) + ((\mathbf{u} \cdot \nabla)A, \psi) - (\nabla \mathbf{u} A, \psi) - (A, \psi \nabla \mathbf{u}) \\ & \quad = -\frac{1}{We}(A, \psi) + \frac{1}{We}(I, \psi). \end{aligned} \quad (2.2)$$

We take $\mathbf{v} = \mathbf{u}$, $q = p$, and keeping in mind that $((\mathbf{u} \cdot \nabla)\mathbf{u}, \mathbf{u}) = 0$, we obtain from the first two equations of System (2.2):

$$\begin{aligned} & \frac{\mathrm{Re}}{2} \frac{d}{dt} \|\mathbf{u}\|^2 + 2\beta \|D\mathbf{u}\|^2 + \frac{1-\beta}{We} \int_{\Omega} A : D\mathbf{u} \\ & \quad + \int_{\Gamma} \alpha |\mathbf{u}_{\tau}|^2 = \int_{\Gamma} \alpha \mathbf{u}_{\tau} \cdot \mathbf{w}_{\tau} + \int_{\Omega} \mathbf{b} \cdot \mathbf{u}. \end{aligned} \quad (2.3)$$
We take $\psi = I$ in $(2.2)_3$ to get:

$$ \frac{d}{dt} \int_{\Omega} \mathrm{tr} \, A + \frac{1}{We} \int_{\Omega} \mathrm{tr} \, A + \int_{\Omega} (\mathbf{u} \cdot \nabla) \mathrm{tr} \, A \\ = 2 \int_{\Omega} A : D\mathbf{u} + \frac{1}{We} d|\Omega|. $$

We use the following Green's formula

$$ \int_{\Omega} (\mathbf{u} \cdot \nabla) \mathrm{tr} \, A = - \int_{\Omega} \mathrm{div} \, \mathbf{u} \mathrm{tr} \, A + \int_{\Gamma} \mathbf{u} \cdot n \mathrm{tr} \, A = 0, $$

and deduce that

$$ \frac{1 - \beta}{2We} \frac{d}{dt} \int \mathrm{tr} \, A + \frac{1 - \beta}{2We^2} \int \mathrm{tr} \, A - \frac{1 - \beta}{2We^2} d|\Omega| \\ = \frac{1 - \beta}{We} \int_{\Omega} A : D\mathbf{u}. \qquad (2.4) $$

We recall that the Korn's inequality reads [10, Chap 3]: there exists $c$ such that for all $v \in \mathbf{V}_{\tau}$,

$$ \int_{\Omega} |Dv|^2 dx + \int_{\Gamma} |v_{\tau}|^2 d\sigma \geq c \int_{\Omega} |\nabla v|^2 dx. \quad (2.5) $$

We replace (2.4) in (2.3), apply Poincaré-Fredrichs's inequality, (2.5), Young's inequality, to obtain

$$ \begin{aligned} & \frac{1}{2} \frac{d}{dt} \left[ \mathrm{Re} \| \mathbf{u} \|^2 + \frac{1-\beta}{We} \int_{\Omega} \mathrm{tr} \, A \right] + 2\beta \| D\mathbf{u} \|^2 \\ & \qquad + \frac{1-\beta}{2We^2} \int_{\Omega} \mathrm{tr} \, A + \int_{\Gamma} \alpha |\mathbf{u}_{\tau}|^2 \\ & = \int_{\Gamma} \alpha^{1/2} \mathbf{u}_{\tau} \cdot \alpha^{1/2} \mathbf{w}_{\tau} + \int_{\Omega} \mathbf{b} \cdot \mathbf{u} + \frac{1-\beta}{2We^2} d|\Omega| \\ & \leq \frac{1}{2} \int_{\Gamma} \alpha |\mathbf{u}_{\tau}|^2 + \frac{1}{2} \int_{\Gamma} \alpha |\mathbf{w}_{\tau}|^2 \\ & \qquad + c \| \mathbf{b} \| (\| \mathbf{u}_{\tau} \|_{\Gamma}^2 + \| D\mathbf{u} \|^2)^{1/2} + \frac{1-\beta}{2We^2} d|\Omega| \\ & \leq \frac{1}{2} \int_{\Gamma} \alpha |\mathbf{u}_{\tau}|^2 + \frac{1}{2} \int_{\Gamma} \alpha |\mathbf{w}_{\tau}|^2 \\ & \qquad + \frac{c}{2\beta} \| \mathbf{b} \|^2 + \frac{\beta}{2} \| \mathbf{u}_{\tau} \|_{\Gamma}^2 + \frac{\beta}{2} \| D\mathbf{u} \|^2 + \frac{1-\beta}{2We^2} d|\Omega|, \end{aligned} $$

which is re-written as follows

$$ \begin{aligned} & \frac{1}{2} \frac{d}{dt} \left[ \mathrm{Re} \| \mathbf{u} \|^2 + \frac{1-\beta}{We} \int_{\Omega} \mathrm{tr} \, A \right] + \frac{3\beta}{2} \| D\mathbf{u} \|^2 \\ & \quad + \frac{1}{2} (\delta - \beta) \| \mathbf{u}_{\tau} \|_{\Gamma}^2 + \frac{1-\beta}{2We^2} \int_{\Omega} \mathrm{tr} \, A \\ & \leq \frac{1}{2} \int_{\Gamma} \alpha |\mathbf{w}_{\tau}|^2 + \frac{c}{2\beta} \| \mathbf{b} \|^2 + \frac{1-\beta}{2We^2} d|\Omega|. \end{aligned} \qquad (2.6) $$

Thus taking $\beta < \delta$, applying (2.5) and Poincare-Fredrichs's inequality again, (2.6) gives

$$ \begin{aligned} & \frac{d}{dt} \left[ \mathrm{Re} \| \mathbf{u} \|^2 + \frac{1-\beta}{We} \int_{\Omega} \mathrm{tr} \, A \right] + c \| \mathbf{u} \|_{\Gamma}^2 \\ & \quad + \frac{1-\beta}{We^2} \int_{\Omega} \mathrm{tr} \, A \\ & \leq \int_{\Gamma} \alpha |\mathbf{w}_{\tau}|^2 + \frac{c}{\beta} \| \mathbf{b} \|^2 + \frac{1-\beta}{We^2} d|\Omega|. \end{aligned} $$

We use the Lemma 2.2 and get:

$$ \begin{aligned} & \frac{d}{dt} \left[ \mathrm{Re} \| \mathbf{u} \|^2 + \frac{1-\beta}{We} \int_{\Omega} \mathrm{tr} \, A \right] \\ & \quad + \min \left( \frac{c}{Re}, \frac{1}{We} \right) \left[ \mathrm{Re} \| \mathbf{u} \|^2 + \frac{1-\beta}{We} \int_{\Omega} \mathrm{tr} \, A \right] \\ & \leq \int_{\Gamma} \alpha |\mathbf{w}_{\tau}|^2 + \frac{c}{\beta} \| \mathbf{b} \|^2 + \frac{1-\beta}{We^2} d|\Omega|. \end{aligned} \qquad (2.7) $$

The result is obtained by using Lemma 2.2 in (2.7). $\square$

**Remark 2.1.** Owing to the relation $E = \frac{1-\beta}{We}(-I + A)$, one deduces that

$$ \begin{aligned} & \mathrm{Re} \| \mathbf{u} \|^2 + \int_{\Omega} \mathrm{tr} \, E + d|\Omega| \frac{1-\beta}{We} \\ & \leq \exp(-\gamma t) \left[ \mathrm{Re} \| \mathbf{u}_0 \|^2 + \int_{\Omega} \mathrm{tr} \, E_0 + d|\Omega| \frac{1-\beta}{We} \right] \\ & \quad + \frac{1-\beta}{We^2} d|\Omega| \frac{1}{\gamma} (1 - \exp(-\gamma t)) \\ & \quad + \frac{1}{\gamma} \| \alpha \|_{L^{\infty}(\Gamma)} \| \mathbf{w}_{\tau} \|_{L^{\infty}(0,t;L^2(\Gamma))}^2 (1 - \exp(-\gamma t)) \\ & \quad + \frac{c}{\gamma \beta} (1 - \exp(-\gamma t)) \| \mathbf{b} \|_{L^{\infty}(0,t;L^2(\Omega))}^2. \end{aligned} $$

**Remark 2.2.** It is apparent that $\mathrm{Re} \| \mathbf{u} \|^2 + \int_{\Omega} \mathrm{tr} \, E dx$ decays slowly when $\mathrm{Re}$ or $\mathrm{We}$ is bigger.

## 3 Space and Time Discretizations

In this section, we propose a discrete scheme associated to the problem described by the equations (1.1)$_1$, (1.1)$_2$, (1.2),(1.4) and (2.1). The aim is to design numerical scheme capable of replicating the energy property display in Proposition 2.1. This will be achieved by adopting the following strategy:

(i) We decompose the conformation tensor $A$ (because it is symmetric and positive definite) and re-write the transport equation,

(ii) we perform space-time approximation.
At the discrete level and for the simplicity, we limit our study to the case $d = 2$. We use the finite element method for the space discretization and the Euler method for the time discretization.

## 3.1 Finite element approximation

In order to use a finite element approximation capable of reproducing the discrete version of Proposition 2.1, we assume that $\Omega$ is a polygon, so it can be completely meshed. We assume also that the triangulation is regular in the sense defined in ([11]): the mesh $(\mathcal{T}_h)_h$ of $\Omega$ is a set of closed non-degenerate triangles called elements, satisfying,

$$ (a) \quad \bar{\Omega} = \bigcup_{1 \le n \le N} K_n . $$

(b) the intersection of two different elements is either empty, a corner, or a whole edge of both elements.

(c) The ratio of the diameter $h_K$ of an element $K$ in $\mathcal{T}_h$ to the diameter of its inscribed circle or sphere is bounded by a constant independent of $K$ and $h$.

As usual, $h$ stands for the maximal diameter of all elements of $\mathcal{T}_h$. For each non-negative integer $l$ and any $K \in \mathcal{T}_h$, $\mathbb{P}_l(K)$ denotes the space of polynomials in $d$ variables restricted to $K$ with total degree less than or equal to $l$.

We consider the conforming approximation of $\mathbb{V}_\tau, M$ and $\mathbb{Z}$ defined as follows:

$$ \mathbb{V}_{\tau,h} = \{ \mathbf{v}_h \in \mathbb{V} \cap \mathcal{C}^0(\bar{\Omega})^2, \text{ for all } K \in \mathcal{T}_{2h}, \\ \mathbf{v}_h|_K \in \mathcal{P}_2(K)^2 \}, $$

$$ L_h^2 = \{ q_h | q_h \in \mathcal{C}^0(\bar{\Omega}), q_h|_K \in \mathcal{P}_1(K), \\ \text{ for all } K \in \mathcal{T}_h \}, $$

$$ M_h = \{ q_h | q_h \in L_h^2, \int_{\Omega} q_h dx = 0 \}, $$

$$ \mathbb{Z}_h = \{ \mathbf{A}_h | \mathbf{A}_h = \begin{pmatrix} a_{1,h} & a_{2,h} \\ a_{2,h} & a_{3,h} \end{pmatrix}, \\ a_i|_K \in \mathcal{P}_1(K), i = 1, 2, 3, \text{ for all } K \in \mathcal{T}_h \}. $$

To simplify the systems of equations established later, we define the following space:

$$ \mathbb{W}_h = \mathbb{V}_{\tau,h} \times M_h \times \mathbb{Z}_h. $$

This is the well-known conforming Taylor-Hood element of degree two with continuous pressures. It should be noted that the degrees of freedom for the velocity are located at the vertices and midpoints of each $K \in \mathcal{T}_h$, while the degrees of freedom for the pressures are associated with the vertices of each triangle $K \in \mathcal{T}_h$.

**Remark 3.1.** It is important to note that there are many possible choices for finite element approximations of the velocity, pressure, and conformation tensor. Readers are invited to consult [8, 31, 32] for further details.

We should also bear in mind that the compatibility condition between $\mathbb{V}_{\tau,h}$ and $M_h$ must hold; that is, there exists $\alpha > 0$ (independent of $h$) such that for all $q_h \in M_h$,

$$ \sup_{0 \neq \mathbf{v}_h \in \mathbb{V}_{\tau,h}} \frac{(q_h, \operatorname{div} \mathbf{v}_h)}{\|\mathbf{v}_h\|_1} \geq \alpha \|q_h\|. \quad (3.1) $$

With the above spaces it is quite natural to approximate problem (2.2) (with obvious notation) by: For a.e $t > 0$, find $(\mathbf{u}_h(t), p_h(t), \mathbf{A}_h(t)) \in \mathbb{W}_h$ such that for all $(\mathbf{v}_h, q_h, \psi_h) \in \mathbb{W}_h$

$$ \begin{aligned} & \mathrm{Re} (\partial_t \mathbf{u}_h, \mathbf{v}_h) + \mathrm{Re} d(\mathbf{u}_h, \mathbf{u}_h, \mathbf{v}_h) + 2\beta(D\mathbf{u}_h, D\mathbf{v}_h) \\ & - (p_h, \operatorname{div} \mathbf{v}_h) + \frac{1-\beta}{We} (\mathbf{A}_h, \nabla \mathbf{v}_h) \\ & = -(\alpha(\mathbf{u}_{\tau,h} - \mathbf{w}_{\tau,h}), \mathbf{v}_{\tau,h})_{\Gamma} + (\mathbf{b}_h, \mathbf{v}_h) \end{aligned} $$

$$ (q_h, \operatorname{div} \mathbf{u}_h) = 0, $$

$$ \begin{aligned} & (\partial_t \mathbf{A}_h, \psi_h) + \tilde{d}(\mathbf{u}_h, \mathbf{A}_h, \psi_h) - (\nabla \mathbf{u}_h \mathbf{A}_h, \psi_h) \\ & - (\mathbf{A}_h, \psi_h \nabla \mathbf{u}_h) = -\frac{1}{We} (\mathbf{A}_h, \psi_h) + \frac{1}{We} (\mathbf{I}, \psi_h), \end{aligned} $$

$$ \mathbf{u}_h(\mathbf{x}, 0) = \mathbf{u}_{0h}, \quad \mathbf{A}_h(\mathbf{x}, 0) = \frac{We}{1-\beta} \mathbf{E}_{0h} + \mathbf{I}, \quad (3.2) $$

with $\mathbf{w}_{\tau,h}$ being an approximation of $\mathbf{w}_{\tau}$ satisfying $\mathbf{w}_{\tau,h} \cdot \mathbf{n}|_{\Gamma} = 0$, $\mathbf{b}_h$ an approximation of $\mathbf{b}$, and $\mathbf{u}_{0h}$ the approximation of $\mathbf{u}_0$ such that $(q_h, \operatorname{div} \mathbf{u}_{0h}) = 0$ for all $q_h \in M_h$. The trilinear forms $d(\cdot, \cdot, \cdot)$ and $\tilde{d}(\cdot, \cdot, \cdot)$ are defined by

$$ \begin{aligned} d(\mathbf{u}, \mathbf{v}, \mathbf{w}) &= \int_{\Omega} (\mathbf{u} \cdot \nabla) \mathbf{v} \cdot \mathbf{w} + \frac{1}{2} \int_{\Omega} \operatorname{div} \mathbf{u} \mathbf{v} \cdot \mathbf{w} \\ &= \frac{1}{2} \left( \int_{\Omega} (\mathbf{u} \cdot \nabla) \mathbf{v} \cdot \mathbf{w} - \int_{\Omega} (\mathbf{u} \cdot \nabla) \mathbf{w} \cdot \mathbf{v} \right), \\ \tilde{d}(\mathbf{u}, \phi, \psi) &= \int_{\Omega} (\mathbf{u} \cdot \nabla) \phi : \psi + \frac{1}{2} \int_{\Omega} \operatorname{div} \mathbf{u} \phi : \psi \\ &= \frac{1}{2} \left( \int_{\Omega} (\mathbf{u} \cdot \nabla) \phi : \psi - \int_{\Omega} (\mathbf{u} \cdot \nabla) \psi : \phi \right), \end{aligned} $$

which display the skew-symmetry property $d(\mathbf{u}, \mathbf{v}, \mathbf{v}) = \tilde{d}(\mathbf{u}, \psi, \psi) = 0$, even when $\mathbf{u}$ does not satisfy the incompressibility condition pointwise but satisfies $\mathbf{u} \cdot \mathbf{n}|_{\partial\Omega} = 0$. We note that
$\tilde{d}(\cdot, \cdot, \cdot)$ and $d(\cdot, \cdot, \cdot)$ are consistent in the sense that

for all $(\phi, \psi) \in H^1(\Omega)^{2 \times 2} \times H^1(\Omega)^{2 \times 2}$ and
$(u, v, w) \in H^1(\Omega) \cap \{\operatorname{div} u = 0\} \times H^1(\Omega)^2$

$$d(u, v, w) = \int_{\Omega} (u \cdot \nabla) v \cdot w dx,$$

$$\tilde{d}(u, \phi, \psi) = \int_{\Omega} (u \cdot \nabla) \phi : \psi dx.$$

It is noted that $\tilde{d}(u_h, A_h, \psi) \equiv 0$ because the extra stress is approximated by piecewise constant function in each element. It is observed in [19] that the backward Euler scheme associated with (3.2) does not replicate the property displayed in Proposition 2.1. In this paper, we are also interested in numerical schemes that preserve the symmetry and positive definiteness of the conformation tensor while ensuring that the discrete level also exhibits the property described in Proposition 2.1. The starting observation is that if the conformation tensor is symmetric and positive definite, then this requirement should be incorporated into the formulation. However, at the discrete level, this constraint is not necessarily satisfied unless it is explicitly enforced by the numerical method. We discuss in the next paragraph how this requirement can be integrated into the formulation.

## 3.2 A reformulation of problem (3.2)

This paragraph follows from [19] and is included here for clarity and a better understanding of what follows. The starting point of this analysis is a simple argument in linear algebra. To reformulate the constitutive equations satisfied by the extra-stress tensor, we observe that $A$ is symmetric and positive definite. From the Cholesky decomposition, it follows that there exists a lower triangular matrix $B$ such that

$$A = BB^T, \qquad (3.3)$$

$B^T$ the transpose of $B$. We replace (2.1) by

$$B_t + (u \cdot \nabla)B - \nabla uB + \frac{1}{2We}B = \frac{1}{2We}B^{-T}, \quad (3.4)$$

and recover $A$ with (3.3). Indeed, we claim that

**Lemma 3.1.** If $B$ solves (3.4), then $A$ solves (2.1).

**Proof.** We take the transpose of (3.4) and multiply by $B$ to obtain

$$BB_t^T + B(u \cdot \nabla)B^T - BB^T(\nabla u)^T + \frac{1}{2We}BB^T = \frac{1}{2We}I. \quad (3.5)$$

Now, (3.5) and (3.4) give

$$(BB^T)_t + (u \cdot \nabla)(BB^T) - \nabla u(BB^T) - BB^T(\nabla u)^T + \frac{1}{We}BB^T = \frac{1}{We}I,$$

which can be written as following:

$$A_t + (u \cdot \nabla)A - \nabla uA - A(\nabla u)^T + \frac{1}{We}A = \frac{1}{We}I.$$

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

Having Lemma 3.1 in mind, we consider the following finite element approximations: For a.e $t > 0$, find $(u_h(t), p_h(t), B_h(t), A_h(t)) \in V_{\tau,h} \times M_h \times X_h \times \mathbb{Z}_h$ such that for all $(v_h, q_h, \psi_h) \in W_h$,

$$\begin{aligned} & \mathrm{Re}(\partial_t u_h, v_h) + \mathrm{Red}(u_h, u_h, v_h) + 2\beta(Du_h, Dv_h) \\ & - (p_h, \operatorname{div} v_h) + \frac{1-\beta}{We}(B_h B_h^T, \nabla v_h) \\ & = -(\alpha(u_{\tau,h} - w_{\tau}), v_{\tau,h})_\Gamma + (b_h, v_h), \end{aligned}$$

$$(q_h, \operatorname{div} u_h) = 0,$$

$$(\partial_t B_h, \psi_h) + \tilde{d}(u_h, B_h, \psi_h) - (\nabla u_h B_h, \psi_h) + \frac{1}{2We}(B_h, \psi_h) - \frac{1}{2We}(B_h^{-T}, \psi_h) = 0$$

$$\begin{aligned} & u_h(x, 0) = u_0, \quad A_h = B_h B_h^T \\ & \text{and} \\ & B_h(x, 0)B_h^T(x, 0) = \frac{We}{1-\beta}E_0 + I, \end{aligned} \quad (3.6)$$

with

$$X_h = \{ B_h | B_h = \begin{pmatrix} a_{1,h} & 0 \\ a_{2,h} & a_{3,h} \end{pmatrix}, \\ a_i | K \in \mathcal{P}_1(K), i = 1, 2, 3, \text{ for all } K \in \mathcal{T}_h \}.$$

We note that (3.6) is nonlinear, and the presence of $B_h^{-T}$ makes it even more complicated and expensive to solve using a direct Euler backward/forward scheme. To overcome this difficulty, and keeping (3.3) in mind, we next formulate a time discretization of problem (3.6) using a semi-implicit scheme.

## 3.3 Time discretization of problem (3.6): Semi-implicit method

In this paragraph, we formulate a semi-implicit scheme based on a time approximation of (3.6) that
preserves the property stated in Proposition 2.1. Then, we study the corresponding properties namely the replication of the energy property.

For the time discretization, we denote by $k$ the difference between two consecutive time points, $t_m$ and $t_{m+1}$, and set $t_m = mk$ for $m \ge 1$. $(\mathbf{u}_{k,h}^m, p_{k,h}^m, A_{k,h}^m, B_{k,h}^m)$ represent the approximations of $(\mathbf{u}_h(t_m), p_h(t_m), A_h(t_m), B_h(t_m))$ at $t_m$.

We note that in [19], the authors introduce a fully implicit scheme associated to (3.6) and given as follows:

**Initialization:** Let $(\mathbf{u}_{k,h}^0, A_{k,h}^0)$ be given, and compute $B_{k,h}^0$ such that $A_{k,h}^0 = B_{k,h}^0 B_{k,h}^{0,T}$.

**Step**($m+1$): Given $(\mathbf{u}_{k,h}^m, A_{k,h}^m, B_{k,h}^m)$, find $(\mathbf{u}_{k,h}^{m+1}, p_{k,h}^{m+1}, B_{k,h}^{m+1}) \in \mathbb{W}_h$

such that for all $(\mathbf{v}_h, q_h, \psi_h) \in \mathbb{W}_h$

$$ \frac{\mathrm{Re}}{k} (\mathbf{u}_{k,h}^{m+1} - \mathbf{u}_{k,h}^m, \mathbf{v}_h) + \mathrm{Red}(\mathbf{u}_{k,h}^{m+1}, \mathbf{u}_{k,h}^{m+1}, \mathbf{v}_h) + 2\beta (D\mathbf{u}_{k,h}^{m+1}, D\mathbf{v}_h) - (\mathbf{p}_{k,h}^{m+1}, \mathrm{div}\,\mathbf{v}_h) + \frac{1-\beta}{We} (\mathbf{B}_{k,h}^{m+1}\mathbf{B}_{k,h}^{m+1,T}, \nabla\mathbf{v}_h) + (\alpha\mathbf{u}_{\tau,k,h}^{m+1}, \mathbf{v}_{\tau,h})_{\Gamma} = (\mathbf{b}_{k,h}^{m+1}, \mathbf{v}_h) + (\alpha\mathbf{w}_{\tau}, \mathbf{v}_{\tau,h})_{\Gamma}, $$

$$ (q_h, \mathrm{div}\,\mathbf{u}_{k,h}^{m+1}) = 0, $$

$$ \frac{1}{k} (\mathbf{B}_{k,h}^{m+1}, \psi_h) + \tilde{d}(\mathbf{u}_{k,h}^{m+1}, \mathbf{B}_{k,h}^{m+1}, \psi_h) - (\nabla \mathbf{u}_{k,h}^{m+1} \mathbf{B}_{k,h}^{m+1}, \psi_h) + \frac{1}{2We} (\mathbf{B}_{k,h}^{m+1}, \psi_h) = \frac{1}{k} (\sqrt{\mathbf{A}_{k,h}^m}, \psi_h), $$

$$ \mathbf{A}_{k,h}^{m+1} = \mathbf{B}_{k,h}^{m+1} \mathbf{B}_{k,h}^{m+1,T} + \frac{k}{We} \mathbf{I}. \qquad (3.7) $$

It is shown in [19] that the scheme (3.7) replicates the energy property. We note that the scheme (3.7) is fully implicit and needs to be linearized at each time step, which means solving an additional inertive system at each time iteration. To avoid this extra iterative system, a new semi-implicit scheme associated with System (3.6) will be introduced and studied in this work.

First we recall the discrete counterpart of Lemma 2.2 (see [36]), which is stated as follows:

**Lemma 3.2.** Let $K > 0$ and let $(\alpha_n)_n, (b_n)_n$ be two sequences of non-negative real numbers satisfying

$$ \forall n, \quad \alpha_n \le (1 + K)\alpha_{n-1} + b_n. $$

Then for all $n \ge 1$, $\alpha_n$ satisfies the bound

$$ \alpha_n \le (1 + K)^n \alpha_0 + \sum_{i=1}^{n} b_i (1 + K)^{n-i}. $$

Next, using the finite element space $\mathbb{W}_h = \mathbb{W}_{\tau,h} \times M_h \times \mathbb{X}_h$, we approximate System (3.6) with the following explicit discrete scheme (using obvious notation):

**Initialization:** Let $(\mathbf{u}_{k,h}^0, A_{k,h}^0)$ be given, and compute $B_{k,h}^0$ such that $A_{k,h}^0 = B_{k,h}^0 B_{k,h}^{0,T}$.

**Step**($m+1$): Given $(\mathbf{u}_{k,h}^m, A_{k,h}^m, B_{k,h}^m)$, find $(\mathbf{u}_{k,h}^{m+1}, p_{k,h}^{m+1}, B_{k,h}^{m+1}) \in \mathbb{W}_h$ such that for all $(\mathbf{v}_h, q_h, \psi_h) \in \mathbb{W}_h$

$$ \frac{\mathrm{Re}}{k} (\mathbf{u}_{k,h}^{m+1} - \mathbf{u}_{k,h}^m, \mathbf{v}_h) + \mathrm{Red}(\mathbf{u}_{k,h}^m, \mathbf{u}_{k,h}^{m+1}, \mathbf{v}_h) + 2\beta (D\mathbf{u}_{k,h}^{m+1}, D\mathbf{v}_h) - (\mathbf{p}_{k,h}^{m+1}, \mathrm{div}\,\mathbf{v}_h) + \frac{1-\beta}{We} (\mathbf{B}_{k,h}^{m+1}\mathbf{B}_{k,h}^{m+1,T}, \nabla\mathbf{v}_h) + (\alpha\mathbf{u}_{\tau,k,h}^{m+1}, \mathbf{v}_{\tau,h})_{\Gamma} = (\mathbf{b}_{k,h}^{m+1}, \mathbf{v}_h) + (\alpha\mathbf{w}_{\tau}, \mathbf{v}_{\tau,h})_{\Gamma}, \\ (q_h, \mathrm{div}\,\mathbf{u}_{k,h}^{m+1}) = 0, \\ \frac{1}{k} (\mathbf{B}_{k,h}^{m+1}, \psi_h) + \tilde{d}(\mathbf{u}_{k,h}^{m+1}, \mathbf{B}_{k,h}^{m+1}, \psi_h) - (\nabla \mathbf{u}_{k,h}^{m+1} \mathbf{B}_{k,h}^{m+1}, \psi_h) + \frac{1}{2We} (\mathbf{B}_{k,h}^{m+1}, \psi_h) = \frac{1}{k} (\sqrt{\mathbf{A}_{k,h}^m}, \psi_h), \\ \mathbf{A}_{k,h}^{m+1} = \mathbf{B}_{k,h}^{m+1} \mathbf{B}_{k,h}^{m+1,T} + \frac{k}{We} \mathbf{I}. \qquad (3.8) $$

Our aim is to show that the solution of the last discrete system satisfies the counterpart of the a priori estimate stated in Proposition 2.1.

The next result establishes the consistency of the semi-implicit scheme (3.8) following similar result in [19].

**Lemma 3.3.** Let $(\mathbf{u}_h, p_h, A_h)$ be solution of (3.2). Then the numerical scheme (3.8) is consistent with (3.2) when $k \to 0$.
**Proof.** With the standard notation, $\mathbf{u}^n \approx \mathbf{u}(t_n)$, we have from the fourth equation of (3.8) that

$$A(t_{m+1}) \approx B(t_{m+1})B^T(t_{m+1}) \quad \text{when} \quad k \to 0. \qquad (3.9)$$

Using Taylor's expansion (having in mind $k \to 0$) one gets

$$B(t_{m+1})B^T(t_m) \approx B(t_{m+1})B^T(t_{m+1}) \approx A(t_{m+1}), \\ d(\mathbf{u}(t_m), \mathbf{u}(t_{m+1}), \mathbf{v}) \approx d(\mathbf{u}(t_{m+1}), \mathbf{u}(t_{m+1}), \mathbf{v}). \qquad (3.10)$$

From (3.9), and (3.10), we deduce that the first equation of (3.8) is consistent with the first equation of (3.2). The third equation of (3.8) is re-written as follows

$$\frac{1}{k} B_{k,h}^{m+1} + \tilde{D}(\mathbf{u}_{k,h}^m, B_{k,h}^{m+1}) - \nabla \mathbf{u}_{k,h}^{m+1} B_{k,h}^m + \frac{1}{2We} B_{k,h}^{m+1} = \frac{1}{k} \sqrt{A_{k,h}^m}, \qquad (3.11)$$

with $\tilde{D}$ given by

$$\tilde{D}(\mathbf{u}_{k,h}^m, B_{k,h}^{m+1}) \phi = \tilde{d}(\mathbf{u}_{k,h}^m, B_{k,h}^{m+1}, \phi).$$

We take the transpose of (3.11) and obtain that

$$\frac{1}{k} B_{k,h}^{m+1,T} + \tilde{D}^T(\mathbf{u}_{k,h}^m, B_{k,h}^{m+1}) - B_{k,h}^{m,T}(\nabla \mathbf{u}_{k,h}^{m+1})^T + \frac{1}{2We} B_{k,h}^{m+1,T} = \frac{1}{k} \left(\sqrt{A_{k,h}^m}\right)^T.$$

At this stage, we continue as in [19, Lemma 1]. $\square$

**Remark 3.2.**

(a) We start the procedure (3.8) with $A_{k,h}^0$ symmetric and positive definite so that $\sqrt{A_{k,h}^0}$ is well defined.

(b) We note from (3.8)₄ that $A_{k,h}^m$ is symmetric and positive definite for all $k$ and $m \ge 0$.

Since the scheme (3.8) is an algebraic system of linear equations, existence and uniqueness of solution are equivalent. We claim that

**Proposition 3.1.** The system of equations (3.8) is feasible.

**Proof.** It suffice to show the uniqueness of the solution of System (3.8), which equivalent to show

that zero is the unique solution of the system:

$$\frac{\mathrm{Re}}{k} (\mathbf{u}_{k,h}^{m+1}, \mathbf{v}_h) + \mathrm{Re} d(\mathbf{u}_{k,h}^m, \mathbf{u}_{k,h}^{m+1}, \mathbf{v}_h) + 2\beta (D\mathbf{u}_{k,h}^{m+1}, D\mathbf{v}_h) - (p_{k,h}^{m+1}, \mathrm{div} \, \mathbf{v}_h) = -\frac{1-\beta}{We} (B_{k,h}^{m+1} B_{k,h}^{m,T}, \nabla \mathbf{v}_h) - (\alpha \mathbf{u}_{\tau,k,h}^{m+1}, \mathbf{v}_{\tau,h})_{\Gamma}, \\ (q_h, \mathrm{div} \, \mathbf{u}_{k,h}^{m+1}) = 0,$$

$$\frac{1-\beta}{kWe} (B_{k,h}^{m+1}, \psi_h) + \frac{1-\beta}{We} \tilde{d}(\mathbf{u}_{k,h}^m, B_{k,h}^{m+1}, \psi_h) - \frac{1-\beta}{We} (\nabla \mathbf{u}_{k,h}^{m+1} B_{k,h}^m, \psi_h) + \frac{1-\beta}{2We^2} (B_{k,h}^{m+1}, \psi_h) = 0.$$

We take $(\mathbf{v}_h, q_h, \psi_h) = (\mathbf{u}_{k,h}^{m+1}, p_{k,h}^{m+1}, B_{k,h}^{m+1})$, and add the resulting equations. This leads to

$$\frac{\mathrm{Re}}{k} \|\mathbf{u}_{k,h}^{m+1}\|^2 + 2\beta \|D\mathbf{u}_{k,h}^{m+1}\|^2 + \left(\alpha \mathbf{u}_{\tau,k,h}^{m+1}, \mathbf{u}_{\tau,h}^{m+1}\right)_{\Gamma} + \frac{1-\beta}{kWe} \|B_{k,h}^{m+1}\|^2 + \frac{1-\beta}{2We^2} \|B_{k,h}^{m+1}\|^2 = 0.$$

The last equation implies that $\mathbf{u}_{k,h}^{m+1} = \mathbf{0}$ and $B_{k,h}^{m+1} = \mathbf{0}$. Reporting back to the system, one gets $(p_{k,h}^{m+1}, \mathrm{div} \, \mathbf{v}_h) = 0$ for all $\mathbf{v}_h$. But making use of the inf-sup condition (see (3.1)), one has $p_{k,h}^{m+1} = 0$. Thus (3.8) has only one solution for all $k, h$. $\square$

The main purpose of this analysis is to formulate the discrete analog to Proposition 2.1. For that purpose, we state the following proposition.

**Proposition 3.2.** Let $(\mathbf{u}_{k,h}^m, B_{k,h}^m, A_{k,h}^m)$ be the discrete solution of (3.8). Assume that $\alpha$ is bounded from below by $\delta$, and $\beta < \delta < \alpha$. Then there exists a constant $c$ depending on $\Omega$ such that

$$\mathrm{Re} \| \mathbf{u}_{k,h}^m \|^2 + \frac{1-\beta}{We} \int_{\Omega} \mathrm{tr} \, A_{k,h}^m \\ \le \frac{1}{(1+\gamma k)^m} \left[ \mathrm{Re} \| \mathbf{u}_{k,h}^0 \|^2 + \frac{1-\beta}{We} \int_{\Omega} \mathrm{tr} \, A_{k,h}^0 \right] \\ + \frac{1}{\gamma} \left( \| \alpha \|_{L^{\infty}(\Gamma)} \| \mathbf{w}_{\tau,h} \|_{\Gamma}^2 \right. \\ \left. + \frac{c}{\beta} \max_{i \ge 1} \| b_{k,h}^i \|^2 + 2 \left( \frac{1}{k} + \frac{1}{We} \right) \frac{k d(1-\beta)}{We^2} |\Omega| \right).$$
If moreover $k < \gamma^{-1} = \left( \min \left( \frac{c\beta}{Re}, \frac{1}{We} \right) \right)^{-1}$, then

$$
\begin{aligned}
& \mathrm{Re} \|\mathbf{u}_{k,h}^m\|^2 + \frac{1-\beta}{We} \int_{\Omega} \mathrm{tr} \, A_{k,h}^m \\
& \le \exp(-mk\gamma/2) \left[ \mathrm{Re} \|\mathbf{u}_{k,h}^0\|^2 + \frac{1-\beta}{We} \int_{\Omega} \mathrm{tr} \, A_{k,h}^0 \right] \\
& \quad + \frac{1}{\gamma} \left( \|\alpha\|_{L^\infty(\Gamma)} \|\mathbf{w}_{\tau,h}\|_{\Gamma}^2 + \frac{c}{\beta} \max_{i \ge 1} \|\mathbf{b}_{k,h}^i\|^2 \right) \\
& \quad + 2 \left( \frac{1}{k} + \frac{1}{We} \right) \frac{kd(1-\beta)}{We^2} |\Omega|.
\end{aligned}
$$

**Proof.** We take $(v_h, q_h, \psi_h) = (\mathbf{u}_{k,h}^{m+1}, p_{k,h}^{m+1}, \mathbf{B}_{k,h}^{m+1})$ in (3.8). One obtains

$$
\begin{aligned}
& \frac{\mathrm{Re}}{2} \|\mathbf{u}_{k,h}^{m+1}\|^2 - \frac{\mathrm{Re}}{2} \|\mathbf{u}_{k,h}^m\|^2 + \frac{\mathrm{Re}}{2} \|\mathbf{u}_{k,h}^{m+1} - \mathbf{u}_{k,h}^m\|^2 \\
& \quad + k\beta \left\| D\mathbf{u}_{k,h}^{m+1} \right\|^2 + k \left( \alpha \mathbf{u}_{\tau,k,h}^{m+1}, \mathbf{u}_{\tau,k,h}^{m+1} \right)_{\Gamma}
\end{aligned}
$$

$$
= -k \frac{1 - \beta}{We} \left( \mathbf{B}_{k,h}^{m+1} \mathbf{B}_{k,h}^{m,T}, \nabla \mathbf{u}_{k,h}^{m+1} \right) + k \left( \alpha \mathbf{w}_{\tau,h}, \mathbf{u}_{\tau,k,h}^{m+1} \right)_{\Gamma} \\
+ k \left( \mathbf{b}_{k,h}^{m+1}, \mathbf{u}_{k,h}^{m+1} \right),
$$

$$
\frac{1 - \beta}{We} \|\mathbf{B}_{k,h}^{m+1}\|^2 - k \frac{1 - \beta}{We} \left( \nabla \mathbf{u}_{k,h}^{m+1}, \mathbf{B}_{k,h}^{m+1} \mathbf{B}_{k,h}^{m,T} \right) \\
+ \frac{k(1 - \beta)}{2We^2} \|\mathbf{B}_{k,h}^{m+1}\|^2 = \frac{1 - \beta}{We} \left( \sqrt{\mathbf{A}_{k,h}^m}, \mathbf{B}_{k,h}^{m+1} \right),
$$

$$
\mathbf{A}_{k,h}^{m+1} = \mathbf{B}_{k,h}^{m+1} \mathbf{B}_{k,h}^{m+1,T} + \frac{k}{We} \mathbf{I}. \qquad (3.12)
$$

Adding the first two equations of (3.12) (and
maintaining the third equation) gives

$$
\begin{aligned}
& \frac{\mathrm{Re}}{2} \|\mathbf{u}_{k,h}^{m+1}\|^2 + \left(1 + \frac{k}{2We}\right) \frac{1-\beta}{We} \|\mathbf{B}_{k,h}^{m+1}\|^2 \\
& \quad + \frac{\mathrm{Re}}{2} \|\mathbf{u}_{k,h}^{m+1} - \mathbf{u}_{k,h}^m\|^2 \\
& \quad + k\beta \left\| D\mathbf{u}_{k,h}^{m+1} \right\|^2 + k \left( \alpha \mathbf{u}_{\tau,k,h}^{m+1}, \mathbf{u}_{\tau,k,h}^{m+1} \right)_{\Gamma} \\
& = \frac{\mathrm{Re}}{2} \|\mathbf{u}_{k,h}^m\|^2 + k \left( \alpha \mathbf{w}_{\tau,h}, \mathbf{u}_{\tau,k,h}^{m+1} \right)_{\Gamma} + k \left( \mathbf{b}_{k,h}^{m+1}, \mathbf{u}_{k,h}^{m+1} \right) \\
& \quad + \frac{1-\beta}{We} \left( \sqrt{\mathbf{A}_{k,h}^m}, \mathbf{B}_{k,h}^{m+1} \right),
\end{aligned}
$$

$$
\mathbf{A}_{k,h}^{m+1} = \mathbf{B}_{k,h}^{m+1} \mathbf{B}_{k,h}^{m+1,T} + \frac{k}{We} \mathbf{I}. \qquad (3.13)
$$

But, using Cauchy-Shwarz's, Hölder's,
Poincare-Fredrichs's inequality, (2.5), and Young's

inequality yields

$$
\begin{aligned}
& k \left( \alpha \mathbf{w}_{\tau,h}, \mathbf{u}_{\tau,k,h}^{m+1} \right)_{\Gamma} \\
& \leq \frac{k}{2} \int_{\Gamma} \alpha |\mathbf{w}_{\tau,h}|^2 + \frac{k}{2} \int_{\Gamma} \alpha |\mathbf{u}_{\tau,k,h}^{m+1}|^2 , \\
& k \left( \mathbf{b}_{k,h}^{m+1}, \mathbf{u}_{k,h}^{m+1} \right) \\
& \leq k \|\mathbf{b}_{k,h}^{m+1}\| \|\mathbf{u}_{k,h}^{m+1}\| \leq k c \|\mathbf{b}_{k,h}^{m+1}\| \|\nabla \mathbf{u}_{k,h}^{m+1}\| \\
& \leq \frac{kc}{\beta} \|\mathbf{b}_{k,h}^{m+1}\|^2 + \frac{k\beta}{2} \|D\mathbf{u}_{k,h}^{m+1}\|^2 + \frac{k\beta}{2} \|\mathbf{u}_{\tau,k,h}^{m+1}\|_{\Gamma}^2,
\end{aligned}
$$

$$
\begin{aligned}
\left( \sqrt{\mathbf{A}_{k,h}^m}, \mathbf{B}_{k,h}^{m+1} \right) &\le \frac{1}{2} \left\| \sqrt{\mathbf{A}_{k,h}^m} \right\|^2 + \frac{1}{2} \left\| \mathbf{B}_{k,h}^{m+1} \right\|^2 \\
&\le \frac{1}{2} \int_{\Omega} \mathrm{tr} \, \mathbf{A}_{k,h}^m + \frac{1}{2} \left\| \mathbf{B}_{k,h}^{m+1} \right\|^2.
\end{aligned}
$$

Returning to (3.13), one obtains

$$
\begin{aligned}
& + \frac{\mathrm{Re}}{2} \|\mathbf{u}_{k,h}^{m+1} - \mathbf{u}_{k,h}^m\|^2 + \frac{k\beta}{2} \left\| D\mathbf{u}_{k,h}^{m+1} \right\|^2 \\
& + \frac{k}{2} \left( \|\alpha\|_{L^\infty(\Gamma)} - \beta \right) \|\mathbf{u}_{\tau,k,h}^{m+1}\|_{\Gamma}^2 \\
& \leq \frac{\mathrm{Re}}{2} \|\mathbf{u}_{k,h}^m\|^2 + \frac{k}{2} \int_{\Gamma} \alpha |\mathbf{w}_{\tau,h}|^2 + \frac{kc}{\beta} \|\mathbf{b}_{k,h}^{m+1}\|^2 \\
& \quad + \frac{1-\beta}{2We} \int_{\Omega} \mathrm{tr} \, \mathbf{A}_{k,h}^m
\end{aligned}
$$

and

$$
\mathbf{A}_{k,h}^{m+1} = \mathbf{B}_{k,h}^{m+1} \mathbf{B}_{k,h}^{m+1,T} + \frac{k}{We} \mathbf{I}. \qquad (3.14)
$$

We take the trace on the second equation in (3.14) and
obtain that

$$
\mathrm{tr} \, \mathbf{A}_{k,h}^{m+1} = \mathbf{B}_{k,h}^{m+1} : \mathbf{B}_{k,h}^{m+1} + \frac{kd}{We} = \left| \mathbf{B}_{k,h}^{m+1} \right|^2 + \frac{kd}{We}.
$$

We then deduce that

$$
\begin{aligned}
& \frac{\mathrm{Re}}{k} \|\mathbf{u}_{k,h}^{m+1}\|^2 + \left(\frac{1}{k} + \frac{1}{We}\right) \frac{1-\beta}{We} \int_{\Omega} \mathrm{tr} \, \mathbf{A}_{k,h}^{m+1} \\
& \quad + \beta \left\| D\mathbf{u}_{k,h}^{m+1} \right\|^2 + \frac{k}{2}(\delta - \beta) \|\mathbf{u}_{\tau,k,h}^{m+1}\|_{\Gamma}^2 \\
& \leq \frac{\mathrm{Re}}{k} \|\mathbf{u}_{k,h}^{m}\|^2 + \frac{1-\beta}{kWe} \int_{\Omega} \mathrm{tr} \, \mathbf{A}_{k,h}^{m} + \int_{\Gamma} \alpha |\mathbf{w}_{\tau,h}|^2 \\
& \quad + 2\frac{c}{\beta} \|\mathbf{b}_{k,h}^{m+1}\|^2 + \left(\frac{1}{k} + \frac{1}{We}\right) \frac{kd(1-\beta)}{We^2} |\Omega|.
\end{aligned}
\tag{3.15}
$$

By using $\beta < \delta$, relation (2.5) and
Poincare-Fredrich's inequality, (3.15) implies
that

$$
\begin{aligned}
& \left( \mathrm{Re} \| \mathbf{u}_{k,h}^{m+1} \|^2 + \frac{1-\beta}{\mathrm{We}} \int_{\Omega} \mathrm{tr} \, A_{k,h}^{m+1} \right) \\
& \le \left[ \frac{1}{1+\gamma k} \right] \left[ \mathrm{Re} \| \mathbf{u}_{k,h}^m \|^2 + \frac{1-\beta}{\mathrm{We}} \int_{\Omega} \mathrm{tr} \, A_{k,h}^m \right] \\
& \quad + \left[ \frac{k}{1+\gamma k} \right] \left( \|\alpha\|_{L^\infty(\Gamma)} \|\mathbf{w}_{\tau,h}\|_{\Gamma}^2 + \frac{c}{\beta} \|\mathbf{b}_{k,h}^{m+1}\|^2 \right. \\
& \qquad \left. + 2 \left( \frac{1}{k} + \frac{1}{\mathrm{We}} \right) \frac{kd(1-\beta)}{\mathrm{We}^2} |\Omega| \right).
\end{aligned}
\quad (3.16)
$$

Thus by induction (application of Lemma 3.2), (3.16)
gives

$$
\begin{aligned}
& \mathrm{Re} \| \mathbf{u}_{k,h}^m \|^2 + \frac{1-\beta}{\mathrm{We}} \int_{\Omega} \mathrm{tr} \, A_{k,h}^m \\
& \le \frac{1}{(1+\gamma k)^m} \left[ \mathrm{Re} \| \mathbf{u}_{k,h}^0 \|^2 + \frac{1-\beta}{\mathrm{We}} \int_{\Omega} \mathrm{tr} \, A_{k,h}^0 \right] \\
& \quad + k \sum_{i=1}^m \frac{1}{(1+\gamma k)^i} \left( \|\alpha\|_{L^\infty(\Gamma)} \| \mathbf{w}_{\tau,h} \|_{\Gamma}^2 + \frac{c}{\beta} \| \mathbf{b}_{k,h}^i \|^2 \right. \\
& \qquad \left. + 2 \left( \frac{1}{k} + \frac{1}{\mathrm{We}} \right) \frac{kd(1-\beta)}{\mathrm{We}^2} |\Omega| \right) \\
& \le \frac{1}{(1+\gamma k)^m} \left[ \mathrm{Re} \| \mathbf{u}_{k,h}^0 \|^2 + \frac{1-\beta}{\mathrm{We}} \int_{\Omega} \mathrm{tr} \, A_{k,h}^0 \right] \\
& \quad + \frac{1}{\gamma} \left( \|\alpha\|_{L^\infty(\Gamma)} \| \mathbf{w}_{\tau,h} \|_{\Gamma}^2 + \frac{c}{\beta} \max_{i \ge 1} \| \mathbf{b}_{k,h}^i \|^2 \right. \\
& \qquad \left. + 2 \left( \frac{1}{k} + \frac{1}{\mathrm{We}} \right) \frac{kd(1-\beta)}{\mathrm{We}^2} |\Omega| \right).
\end{aligned}
\quad (3.17)
$$

Note that exp(x/2) ≤ x + 1 if 0 < x < 1. Hence for k < 1/γ, one has (1 + γk)^-m < exp(-mkγ/2). Then, (3.17) becomes

$$
\begin{aligned}
& \mathrm{Re} \| \mathbf{u}_{k,h}^m \|^2 + \frac{1-\beta}{\mathrm{We}} \int_{\Omega} \mathrm{tr} \, A_{k,h}^m \\
& \le \exp(-mk\gamma/2) \left[ \mathrm{Re} \| \mathbf{u}_{k,h}^0 \|^2 + \frac{1-\beta}{\mathrm{We}} \int_{\Omega} \mathrm{tr} \, A_{k,h}^0 \right] \\
& \quad + \frac{1}{\gamma} \left( \|\alpha\|_{L^\infty(\Gamma)} \| \mathbf{w}_{\tau,h} \|_{\Gamma}^2 + \frac{c}{\beta} \max_{i \ge 1} \| \mathbf{b}_{k,h}^i \|^2 \right. \\
& \qquad \left. + 2 \left( \frac{1}{k} + \frac{1}{\mathrm{We}} \right) \frac{kd(1-\beta)}{\mathrm{We}^2} |\Omega| \right),
\end{aligned}
$$

and the proof is complete.
□

In the next section, we propose a time discretization
of (3.6) based on operator splitting.

4 Time discretization of problem
(3.6): Marchuk-Yanenko’s scheme

Most current viscoelastic solvers are based on
operator splitting algorithms (see [23, 32, 24]) to
capture specific phenomena at each step and, more
importantly, to simplify computations compared
to the unsplit equations. This approach remains
applicable to the initial value problem (3.6). Indeed,
from an abstract point of view, (3.6) can be rewritten
as follows:

$$
\begin{equation}
\begin{aligned}
& \frac{d\psi}{dt} + A_1(\psi) + A_2(\psi) + A_3(\psi) = \mathbf{f}, \quad (4.1) \\
& \psi(\mathbf{0}) = \psi_0,
\end{aligned}
\end{equation}
$$

where $A_i$ are single-valued operators. To solve
(4.1) numerically, we propose a Marchuk-Yanenko
approach, which is at most first-order accurate.
However, its low accuracy is compensated by good
stability and robustness properties. It is formulated as
follows:

$$
\left\{
\begin{array}{l}
\displaystyle \frac{\psi^{m+1/3} - \psi^m}{\psi^{m+2/3} - \psi^{m+1/3}} + A_1(\psi^{m+1/3}) = \mathbf{f}_1^{m+1}, \\[1em]
\displaystyle \frac{k}{\psi^{m+1} - \psi^{m+2/3}} + A_2(\psi^{m+2/3}) = \mathbf{f}_2^{m+1}, \\[1em]
\displaystyle \frac{\psi^{m+1} - \psi^{m+2/3}}{k} + A_3(\psi^{m+1}) = \mathbf{f}_3^{m+1},
\end{array}
\right.
\qquad (4.2)
$$

with $\sum_{i=1}^3 \mathbf{f}_i^{m+1} = \mathbf{f}^{m+1}$. Applying the scheme (4.2) to the problem (3.6), we obtain:

$$
\mathbf{u}_{k,h}^{0} = \mathbf{u}_{0h}, \quad \mathbf{A}_{k,h}^{0} = \mathbf{A}_{0h},
$$

for $m \ge 0$, we compute $(\mathbf{u}_{k,h}^{m+1/3}, \mathbf{B}_{k,h}^{m+1/3})$,
$(\mathbf{u}_{k,h}^{m+2/3}, \mathbf{p}_{k,h}^{m+2/3}, \mathbf{B}_{k,h}^{m+2/3})$, $(\mathbf{u}_{k,h}^{m+1}, \mathbf{B}_{k,h}^{m+1})$ and
$\mathbf{A}_{k,h}^{m+1}$ via the solution of:

for all $(\psi_h, v_h, q_h) \in \mathbb{Z}_h \times \mathbb{V}_{\tau,h} \times M_h$

Step 1: linear advection equations

$$
\begin{equation}
\begin{split}
& \frac{\mathrm{Re}}{k} (\mathbf{u}_{k,h}^{m+1/3} - \mathbf{u}_{k,h}^{m}, \mathbf{v}_h) \\
& \quad + \mathrm{Re} \, d(\mathbf{u}_{k,h}^{m}, \mathbf{u}_{k,h}^{m+1/3}, \mathbf{v}_h) = 0,
\end{split}
\tag{4.3}
\end{equation}
$$

$$
\frac{1}{k} \left( \mathbf{B}_{k,h}^{m+1/3} - \sqrt{\mathbf{A}_{k,h}^{m}}, \psi_h \right) \\
+ \tilde{d}(\mathbf{u}_{k,h}^{m+1/3}, \mathbf{B}_{k,h}^{m+1/3}, \psi_h) = 0.
$$
**Step 2:** uncoupled linear equations

$$ \frac{Re}{k} (\mathbf{u}_{k,h}^{m+2/3} - \mathbf{u}_{k,h}^{m+1/3}, \mathbf{v}_h) + \beta(D\mathbf{u}_{k,h}^{m+2/3}, D\mathbf{v}_h) - (p_{k,h}^{m+2/3}, \mathrm{div} \, \mathbf{v}_h) = 0, \\ (q_h, \mathrm{div} \, \mathbf{u}_{k,h}^{m+2/3}) = 0, $$

$$ \frac{1}{k} (\mathbf{B}_{k,h}^{m+2/3} - \mathbf{B}_{k,h}^{m+1/3}, \psi_h) + \frac{1}{4We} (\mathbf{B}_{k,h}^{m+2/3}, \psi_h) = 0. \quad (4.4) $$

**Step 3:** coupled linear system

$$ \frac{Re}{k} (\mathbf{u}_{k,h}^{m+1} - \mathbf{u}_{k,h}^{m+2/3}, \mathbf{v}_h) + \beta(D\mathbf{u}_{k,h}^{m+1}, D\mathbf{v}_h) + \frac{1-\beta}{We}(\mathbf{B}_{k,h}^{m+1}\mathbf{B}_{k,h}^{m+2/3,T}, \nabla\mathbf{v}_h) = -(\alpha(\mathbf{u}_{\tau,k,h}^{m+1} - \mathbf{w}_{\tau}), \mathbf{v}_{\tau,h})_{\Gamma} + (\mathbf{b}_{k,h}^{m+1}, \mathbf{v}_h), $$

$$ \frac{1}{k} (\mathbf{B}_{k,h}^{m+1} - \mathbf{B}_{k,h}^{m+2/3}, \psi_h) + \frac{1}{4We} (\mathbf{B}_{k,h}^{m+1}, \psi_h) - (\nabla \mathbf{u}_{k,h}^{m+1} \mathbf{B}_{k,h}^{m+2/3}, \psi_h) = 0, \quad (4.5) $$

and

$$ \mathbf{A}_{k,h}^{m+1} = \mathbf{B}_{k,h}^{m+1} \mathbf{B}_{k,h}^{m+1,T} + \frac{k}{We} \mathbf{I}. \quad (4.6) $$

The scheme (4.3), (4.4), (4.5) and (4.6) allows us to decouple the following difficulties:

1. the advection terms,

2. the incompressibility condition and the related unknown pressure.

**Remark 4.1.** Operator splitting in this context should be viewed as both a time-discretization and an iterative scheme. It is worth noting that other splitting schemes are possible, some of which are certainly more accurate than the one presented above. We believe that the scheme (4.3), (4.4), and (4.5) is one of the simplest schemes involving three operators, designed to satisfy the discrete analog of Proposition 2.1.

**Remark 4.2.** Equation (4.3) is decoupled, and the solution $(\mathbf{u}_{k,h}^{m+1/3}, \mathbf{B}_{k,h}^{m+1/3})$ is well defined. Indeed, these are algebraic linear systems, meaning that existence and uniqueness are equivalent. To prove the uniqueness of the solution $\mathbf{u}_{k,h}^{m+1/3}$, we set $\mathbf{u}_{k,h}^m = \mathbf{0}$

and take $\mathbf{v} = \mathbf{u}_{k,h}^{m+1/3}$, which leads to $\mathbf{u}_{k,h}^{m+1/3} = \mathbf{0}$. Thus $\mathbf{u}_{k,h}^{m+1/3}$ is uniquely defined. Similarly, we can readily check that if $\mathbf{A}_{k,h}^m = \mathbf{0}$, then $\mathbf{B}_{k,h}^{m+1/3} = \mathbf{0}$, ensuring that $\mathbf{B}_{k,h}^{m+1/3} = \mathbf{0}$ is also uniquely defined. The stiffness matrix obtained in this step is not constant.

**Remark 4.3.** Equation (4.4) consists of two completely decoupled systems of equations. First, we solve for $(\mathbf{u}_{k,h}^{m+2/3}, p_{k,h}^{m+2/3})$, which corresponds to a Stokes system for which the existence of a solution is well known (see [16, 7]). The system of equations for $\mathbf{B}_{k,h}^{m+2/3}$ is linear, and its stiffness matrix is positive definite. Hence $\mathbf{B}_{k,h}^{m+2/3}$ is well defined. Moreover, the stiffness matrices obtained for both systems are constant. As a result, they can be factored once, leading to computational time savings.

**Remark 4.4.** Equation (4.5) is linear and coupled. The stiffness matrix is not constant, so the existence of a solution is equivalent to its uniqueness. To establish uniqueness, it suffices to verify that zero is the only solution of

$$ \frac{Re}{k} (\mathbf{u}_{k,h}^{m+1}, \mathbf{v}_h) + \beta(D\mathbf{u}_{k,h}^{m+1}, D\mathbf{v}_h) + \frac{1-\beta}{We}(\mathbf{B}_{k,h}^{m+1}\mathbf{B}_{k,h}^{m+2/3,T}, \nabla\mathbf{v}_h) = -(\alpha\mathbf{u}_{\tau,k,h}^{m+1}, \mathbf{v}_{\tau,h})_{\Gamma}, $$

$$ \frac{1}{k} (\mathbf{B}_{k,h}^{m+1}, \psi_h) + \frac{1}{4We} - (\nabla \mathbf{u}_{k,h}^{m+1} \mathbf{B}_{k,h}^{m+2/3}, \psi_h) = 0. $$

We take $(\mathbf{v}_h, \psi_h) = (\mathbf{u}_{k,h}^{m+1}, \mathbf{B}_{k,h}^{m+1})$, multiply the second equation by $\frac{1-\beta}{We}$ and then add the resulting equations to obtain

$$ \frac{Re}{k} \|\mathbf{u}_{k,h}^{m+1}\|^2 + \beta \|D\mathbf{u}_{k,h}^{m+1}\|^2 + (\alpha \mathbf{u}_{\tau,k,h}^{m+1}, \mathbf{u}_{\tau,k,h}^{m+1})_{\Gamma} + \frac{1-\beta}{kWe} \|\mathbf{B}_{k,h}^{m+1}\|^2 + \frac{1-\beta}{4We^2} \|\mathbf{B}_{k,h}^{m+1}\|^2 = 0. $$

Hence $\|\mathbf{u}_{k,h}^{m+1}\| = \|\mathbf{B}_{k,h}^{m+1}\| = 0$. Thus zero is the solution.

The key property of the scheme (4.3), (4.4), (4.5), and (4.6) is as follows:

**Proposition 4.1.** Let $(\mathbf{u}_{k,h}^m, \mathbf{B}_{k,h}^m, \mathbf{A}_{k,h}^m)$ be given by (4.3), (4.4) and (4.5), (4.6). Assume that $\alpha$ is bounded from below by $\delta$ and $\beta < \delta < \alpha$. Then there exists a
constant c depending on Ω such that

$$
\begin{align*}
& \mathrm{Re} \| \mathbf{u}_{k,h}^m \|^2 + \frac{1 - \beta}{We} \int_{\Omega} \mathrm{tr} \, \mathbf{A}_{k,h}^m \\
& \leq \frac{1}{(1 + \gamma k)^m} \left[ \mathrm{Re} \| \mathbf{u}_{k,h}^0 \|^2 + \frac{1 - \beta}{We} \int_{\Omega} \mathrm{tr} \, \mathbf{A}_{k,h}^0 \right] \\
& \quad + \frac{1}{\gamma} \left( \| \alpha \|_{L^{\infty}(\Gamma)} \| \mathbf{w}_{\tau,h} \|_{\Gamma}^2 + \frac{c}{\beta} \max_{i \geq 1} \| \mathbf{b}_{k,h}^i \|^2 \right. \\
& \qquad \left. + 2 \left( \frac{1}{k} + \frac{1}{We} \right) \frac{kd(1 - \beta)}{We^2} |\Omega| \right).
\end{align*}
$$

If moreover $k < \gamma^{-1} = \left( \min \left( \frac{c\beta}{Re}, \frac{1}{We} \right) \right)^{-1}$, then

$$
\begin{align*}
& \mathrm{Re} \| \mathbf{u}_{k,h}^m \|^2 + \frac{1 - \beta}{We} \int_{\Omega} \mathrm{tr} \, \mathbf{A}_{k,h}^m \\
& \le \exp(-mk\gamma/2) \left[ \mathrm{Re} \| \mathbf{u}_{k,h}^0 \|^2 + \frac{1 - \beta}{We} \int_{\Omega} \mathrm{tr} \, \mathbf{A}_{k,h}^0 \right] \\
& \quad + \frac{1}{\gamma} \left( \|\alpha\|_{L^{\infty}(\Gamma)} \|\mathbf{w}_{\tau,h}\|_{\Gamma}^2 + \frac{c}{\beta} \max_{i \ge 1} \|\mathbf{b}_{k,h}^i\|^2 \right. \\
& \qquad \left. + 2 \left( \frac{1}{k} + \frac{1}{We} \right) \frac{kd(1-\beta)}{We^2} |\Omega| \right).
\end{align*}
$$

**Proof.** We take $(v_h, \psi_h) = 2(\mathbf{u}_{k,h}^{m+1/3}, \mathbf{B}_{k,h}^{m+1/3})$ in (4.3) and obtain

$$
\begin{align*}
& \| \mathbf{u}_{k,h}^{m+1/3} \|^2 + \| \mathbf{u}_{k,h}^{m+1/3} - \mathbf{u}_{k,h}^m \|^2 = \| \mathbf{u}_{k,h}^m \|^2, \\
& \| \mathbf{B}_{k,h}^{m+1/3} \|^2 + \| \mathbf{B}_{k,h}^{m+1/3} - \sqrt{\mathbf{A}_{k,h}^m} \|^2 \\
& \qquad = \| \sqrt{\mathbf{A}_{k,h}^m} \|^2 = \int_{\Omega} \mathrm{tr} \, \mathbf{A}_{k,h}^m,
\end{align*}
$$

from which we deduce that

$$
\begin{align}
\| \mathbf{u}_{k,h}^{m+1/3} \|^2 &\le \| \mathbf{u}_{k,h}^m \|^2, \tag{4.7} \\
\| \mathbf{B}_{k,h}^{m+1/3} \|^2 &\le \int_{\Omega} \mathrm{tr} \, \mathbf{A}_{k,h}^m. \nonumber
\end{align}
$$

We take $(v_h, \psi_h) = 2(\mathbf{u}_{k,h}^{m+2/3}, \mathbf{B}_{k,h}^{m+2/3})$ and $q_h = p_{k,h}^{m+2/3}$ in (4.4), and obtain

$$
\begin{align*}
& \frac{\mathrm{Re}}{k} \left( \| \mathbf{u}_{k,h}^{m+2/3} \|^2 - \| \mathbf{u}_{k,h}^{m+1/3} \|^2 + \| \mathbf{u}_{k,h}^{m+2/3} - \mathbf{u}_{k,h}^{m+1/3} \|^2 \right) \\
& \quad + 2\beta \| D \mathbf{u}_{k,h}^{m+2/3} \|^2 = 0, \\
& \frac{1}{k} \left( \| \mathbf{B}_{k,h}^{m+2/3} \|^2 - \| \mathbf{B}_{k,h}^{m+1/3} \|^2 + \| \mathbf{B}_{k,h}^{m+2/3} - \mathbf{B}_{k,h}^{m+1/3} \|^2 \right) \\
& \quad + \frac{1}{2We} \| \mathbf{B}_{k,h}^{m+2/3} \|^2 = 0,
\end{align*}
$$

from which we deduce that

$$
\begin{align}
& \mathrm{Re} \| \mathbf{u}_{k,h}^{m+2/3} \|^2 + 2k\beta \| D \mathbf{u}_{k,h}^{m+2/3} \|^2 \le \mathrm{Re} \| \mathbf{u}_{k,h}^{m+1/3} \|^2, \nonumber \\
& \qquad \left( 1 + \frac{k}{2We} \right) \| \mathbf{B}_{k,h}^{m+2/3} \|^2 \le \| \mathbf{B}_{k,h}^{m+1/3} \|^2. \tag{4.8}
\end{align}
$$

Replacing (4.7) in (4.8) gives

$$
\begin{align}
& \mathrm{Re} \| \mathbf{u}_{k,h}^{m+2/3} \|^2 + 2k\beta \| D \mathbf{u}_{k,h}^{m+2/3} \|^2 \le \mathrm{Re} \| \mathbf{u}_{k,h}^{m+1/3} \|^2, \nonumber \\
& \frac{1 - \beta}{We} \left( 1 + \frac{k}{2We} \right) \| \mathbf{B}_{k,h}^{m+2/3} \|^2 \le \frac{1 - \beta}{We} \int_{\Omega} \mathrm{tr} \, \mathbf{A}_{k,h}^m. \tag{4.9}
\end{align}
$$

Finally we take $(v_h, \psi_h) = 2(\mathbf{u}_{k,h}^{m+1}, \mathbf{B}_{k,h}^{m+1})$ in (4.5).
One obtains

$$
\begin{align*}
& \mathrm{Re} \left( \| \mathbf{u}_{k,h}^{m+1} \|^2 + \| \mathbf{u}_{k,h}^{m+1} - \mathbf{u}_{k,h}^{m+2/3} \|^2 \right) \\
& \quad + 2k\beta \| D \mathbf{u}_{k,h}^{m+1} \|^2 \\
& \quad + 2k \frac{(1-\beta)}{We} (\mathbf{B}_{k,h}^{m+1} \mathbf{B}_{k,h}^{m+2/3,T}, \nabla \mathbf{u}_{k,h}^{m+1}) \\
& = \mathrm{Re} \| \mathbf{u}_{k,h}^{m+2/3} \|^2 - 2k(\alpha(\mathbf{u}_{k,h,\tau}^{m+1} - \mathbf{w}_{\tau,h}), \mathbf{u}_{\tau,k,h}^{m+1})_{\Gamma} \\
& \quad - 2k(\mathbf{b}_{k,h}^{m+1}, \mathbf{u}_{k,h}^{m+1}), \\
& \frac{1-\beta}{We} \| \mathbf{B}_{k,h}^{m+1} \|^2 + \frac{1-\beta}{We} \| \mathbf{B}_{k,h}^{m+1} - \mathbf{B}_{k,h}^{m+2/3} \|^2 \\
& \quad + \frac{1-\beta}{We} \frac{k}{2We} \| \mathbf{B}_{k,h}^{m+1} \|^2 \\
& \quad - 2k \frac{1-\beta}{We} (\nabla \mathbf{u}_{k,h}^{m+1} \mathbf{B}_{k,h}^{m+2/3}, \mathbf{B}_{k,h}^{m+1}) \\
& = \frac{1-\beta}{We} \| \mathbf{B}_{k,h}^{m+2/3} \|^2.
\end{align*}
$$

Now, we add these relations and discard some positive
terms to get

$$
\begin{align*}
& \mathrm{Re} \| \mathbf{u}_{k,h}^{m+1} \|^2 + 2k\beta \| D \mathbf{u}_{k,h}^{m+1} \|^2 + \frac{1-\beta}{We} \| \mathbf{B}_{k,h}^{m+1} \|^2 \\
& \quad + \frac{k}{2We} \frac{1-\beta}{We} \| \mathbf{B}_{k,h}^{m+1} \|^2 \\
& \quad + 2k \int_{\Gamma} \alpha |\mathbf{u}_{k,h,\tau}^{m+1}|^2 \\
& \leq \mathrm{Re} \| \mathbf{u}_{k,h}^{m+2/3} \|^2 + \frac{1-\beta}{We} \| \mathbf{B}_{k,h}^{m+2/3} \|^2 \\
& \quad + 2k \int_{\Gamma} \alpha \mathbf{u}_{k,h,\tau}^{m+1} \cdot \mathbf{w}_{\tau,h} - 2k(\mathbf{b}_{k,h}^{m+1}, \mathbf{u}_{k,h}^{m+1}) \\
& \leq \mathrm{Re} \| \mathbf{u}_{k,h}^{m+2/3} \|^2 + \frac{1-\beta}{We} \| \mathbf{B}_{k,h}^{m+2/3} \|^2 \\
& \quad + k \int_{\Gamma} \alpha |\mathbf{u}_{k,h,\tau}^{m+1}|^2 + k \int_{\Gamma} \alpha |\mathbf{w}_{\tau,h}|^2 \\
& \quad + 2ck \| \mathbf{b}_{k,h}^{m+1} \| (\| D \mathbf{u}_{k,h}^{m+1} \|^2 + \| \mathbf{u}_{\tau,k,h}^{m+1} \|_{\Gamma}^2)^{1/2},
\end{align*}
$$
which after Young's inequality gives

$$
\begin{align}
& \mathrm{Re} \| \mathbf{u}_{k,h}^{m+1} \|^2 + k\beta \| D \mathbf{u}_{k,h}^{m+1} \|^2 + \frac{1-\beta}{We} \| \mathbf{B}_{k,h}^{m+1} \|^2 \nonumber \\
& + \frac{k}{2We} \frac{1-\beta}{We} \| \mathbf{B}_{k,h}^{m+1} \|^2 + k(\delta - \beta) \| \mathbf{u}_{k,h,\tau}^{m+1} \|^2_{\Gamma} \nonumber \\
& \leq \mathrm{Re} \| \mathbf{u}_{k,h}^{m+2/3} \|^2 + \frac{1-\beta}{We} \| \mathbf{B}_{k,h}^{m+2/3} \|^2 \nonumber \\
& + k \int_{\Gamma} \alpha |\mathbf{w}_{\tau,h}|^2 + \frac{ck}{\beta} \| \mathbf{b}_{k,h}^{m+1} \|^2. \tag{4.10}
\end{align}
$$

Finally, the asserted result is obtained by using (4.9), (4.10), the relation

$$
\int_{\Omega} \mathrm{tr} \, A_{k,h}^{m+1} = \int_{\Omega} |\mathbf{B}_{k,h}^{m+1}|^2 + \frac{kd}{We} |\Omega|,
$$

and lemma 3.2.

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

5 Numerical simulations

In this section, we present numerical simulations
performed using FreeFem++ [25]. The primary
goal is to computationally demonstrate that the
schemes (3.8) and (4.3)-(4.6) preserve the energy
properties stated in Proposition 3.2 and Proposition
4.1, respectively. In the results that follow, the
scheme (3.8) is referred to as the semi-implicit
scheme, while (4.3)-(4.6) is called Scheme 1. We
evaluate the performance of our new schemes using
the lid-driven cavity flow, a well-known benchmark
for flow problems.

We present two test cases. In the first, we consider
a nonzero external source **b**, and in the second, we
examine the lid-driven cavity. In both cases, we
consider the square domain Ω = (0, 1)², where each
edge is divided into N segments of equal length.
Consequently, the corresponding mesh contains 2N²
elements.

**5.1 First test case**

We begin by setting $N = 30$ and considering the
following external force:

$$
\mathbf{b} = \left( \begin{array}{l} (1 + \tanh(t)) x^2 (x-1)^2 y^2 (y-1)^2 \\ e^{-t} xy \end{array} \right).
$$

The initial values of the velocity and the elastic stress are set to **u**₀ = **0** and **E**₀ = 0.2*I*. Additionally, **w**<sub>τ</sub> = **w**·**τ** where **w** = (0.2, 0) on the top and bottom boundaries, and **w** is identically zero on the left and right sides.

We begin by comparing the semi-implicit scheme (3.8) to Scheme 1. Figure 1 illustrates the evolution of the first velocity component $u_1$ at the point (0.25, 0.25) and the energy

$$
E_m = \mathrm{Re} \| \mathbf{u}_{k,h}^m \|^2 + \frac{1 - \beta}{We} \int_{\Omega} \mathrm{tr} \, A_{k,h}^m dx
$$

over time, for $k = 0.002$, $Re = 1$, $We = 1$, $\beta = 0.5$ and $\alpha = 1$, with a final time $T = 100$. The results indicate that the two schemes produce very similar outcomes.

<figure><img src="image_4.png" /><figcaption>Figure 1: Comparison of $u_1$ (up) and $E_m$ (down) between the scheme (3.8) and scheme 1 (k=0.002, $Re=1$, $We=1$, $\beta = 0.5$, $\alpha = 1$)</figcaption></figure>

Next, we examine the dependency of energy with
respect to Re for k = 0.002, We = 1, β =
0.5 and α = 1. Figure 2 shows the evolution
of energy Em with respect to the time for various
values of Re using scheme (3.8). We observe that all
energy curves exhibit a decay and remain similarly
bounded. Moreover, Em decreases more slowly as
the Reynolds number increases, a behavior that is also
supported by Proposition 3.2.

<figure><img src="image_5.png" /><figcaption>Figure 2: $E_m$ for (3.8) with respect to time for $Re=1,5,10,15,20,50,100$ (k=0.002, $We=1$, $\beta = 0.5$, $\alpha = 1$).</figcaption></figure>

Furthermore, to analyze the performance of the
algorithms for large Reynolds numbers, Figure 3
presents a comparison of the energy $E_m$ between (3.8)
and Scheme 1 for Re=1000,2000,3000, N = 60, k = 0.005 and T = 500. We observe that for each Reynolds number, the results obtained with the semi-implicit scheme (3.8) and Scheme 1 remain very close.

<figure><img src="image_6.png" /><figcaption>Figure 3: $E_m$ for (3.8) and scheme 1 with respect to time for Re=1000,2000,3000 (k=0.005, We=1, $\beta$ = 0.5, $\alpha$ = 1).</figcaption></figure>

## 5.2 Lid Driven cavity

In this paragraph, we present the numerical results for the Lid-Driven Cavity Oldroyd-B model under the slip boundary condition (1.1)–(1.5). This is a widely studied example that has been analyzed in various contexts in [30, 31, 32, 1]. The fluid is confined in the $xy$-plane with $(x, y) \in (0, 1)^2$, and the velocity $u$ satisfies $u \cdot n = 0$ on $\Gamma$. On the boundary $\{y = 1, 0 < x < 1\}$ we have

$$ (Tn)_{\tau} = (8(1 + \tanh(t)) x^2 (x - 1)^2, 0)^T \quad (5.1) $$

while the slip boundary conditions $(Tn)_{\tau} + u_{\tau} = 0$ (with $\alpha = 1$ and $w = 0$) is imposed on all other boundaries. Thus, the flow starts smoothly. It is observed that when $t$ is sufficiently large, $(Tn)_{\tau}$ reaches its maximum value $(1, 0)^T$ at the center $x = 1/2$. This boundary condition prevents the formation of local singularities at the top-right and top-left corners. The initial condition for the elastic stress is $E_0 = 0.2t$, and $We=1$. In all the numerical results presented in this section, we set N=60, k=0.005, $\alpha = 1$, $\beta = 0.5$ and $b = 0$.

First, we compare the first element $a_{11}$ of the tensor $A_{k,h}^m$, the first component $u_1$ of the velocity $u = (u_1, u_2)$, and the energy $E_m$ at the point $(0, 25, 0.25)$ between (3.8) and Scheme 1 for $Re = 1$. Figures 4, 5, and 6 demonstrate that in this test, (3.8) and Scheme 1 exhibit the same behavior and are very close.

Furthermore, we show in Figure 7 and 8 the first and second components $(u_1$ and $u_2)$ of the velocity $u$ at the vertical line $x = 0.5$ (with respect to $y$) and for $t = 5$ by using Scheme (3.8).

<figure><img src="image_7.png" /><figcaption>Figure 4: Comparison of $a_{11}$ between the scheme (3.8) and scheme 1 for k=0.005, Re=1, We=1, $\beta$ = 0.5, $\alpha$ = 1</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 5: Comparison of $u_2$ and $E_m$ between the scheme (3.8) and scheme 1 for k=0.005, Re=1, We=1, $\beta$ = 0.5, $\alpha$ = 1</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 6: Comparison of $E_m$ between the scheme (3.8) and scheme 1 for k=0.005, Re=1, We=1, $\beta$ = 0.5, $\alpha$ = 1</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 7: First component $u_1$ of the velocity by using Scheme (3.8) with T=5, k=0.005, Re=1, We=1, $\beta$ = 0.5, $\alpha$ = 1</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 8: Second component $u_2$ of the velocity by using Scheme (3.8) with $T=5$, $k=0.005$, $Re=1$, $We=1$, $\beta = 0.5$, $\alpha = 1$.</figcaption></figure>

Next, we present comparisons for large Reynolds numbers. Figures 9 and 10 show comparisons of the energy $E_m$ between (3.8) and Scheme 1 for $Re=1000, 2000$ and $T=500$. We observe good stability in the algorithms, as the asymptotic results are very close.

<figure><img src="image_12.png" /><figcaption>Figure 9: graph $E_m$ for the scheme (3.8) and scheme 1 with respect to time for $Re=1000$ and $k=0.005$, $We=1$, $\beta = 0.5$, $\alpha = 1$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 10: graph $E_m$ for the scheme (3.8) and scheme 1 with respect to time for $Re=2000$ and $k=0.005$, $We=1$, $\beta = 0.5$, $\alpha = 1$.</figcaption></figure>

Furthermore, Figures 11 show the stream function $\psi \in H_0^1(\Omega)$ such that $u = \text{curl} \psi$ for $Re = 1000$, $t = 5, 10, 20, 50$ with (3.8).

<figure><img src="image_14.png" /><figcaption>Figure 11: stream function $\psi$ for $Re=1000$ and for $t = 1$ (First figure), $t = 10$ (Second figure), $t = 20$ (Third figure) and $t = 50$ (Fourth figure).</figcaption></figure>
**Algorithms comparison** Let us now compare the CPU time between the semi-implicit scheme and Scheme 1. For this, we consider $Re=1,We=1, \beta = 0.5, \alpha = 1, T = 10$ and $k = 0.005$, and compute the CPU time for both algorithms. Table 1 shows that the semi-implicit scheme is faster than Scheme 1. In fact, upon closer inspection of the two algorithms, we observe that the matrices to be assembled are the same at each iteration, but the difference lies in the resolution level. The number of unknowns in (3.8) (semi-implicit scheme) is of the order $15N^2$, while the unknowns for only the third step of Scheme 1 is of the order $14N^2$ (the difference being the pressure). We should note that for Scheme 1, there are also the first two steps where the number of unknowns is relatively large. This is why, in our opinion, in this particular case where the velocity and pressure are approximated with the Taylor-Hood elements, and the stress $P_0$ element, the semi-implicit scheme is faster than Scheme 1.

<table><caption>Table 1: Performances of the Semi-implicit scheme (3.8) and scheme 1 on the Lid Driven for $Re=1,We=1, \beta = 0.5, k = 0.005, \alpha = 1, T = 10$</caption><thead><tr><th rowspan="2">h</th><th>scheme (3.8)</th><th>scheme 1</th><th rowspan="2">The quotient:<br>(Scheme 1)<br>(Semi-implicit)</th></tr><tr><th>CPU time (s)</th><th>CPU time (s)</th></tr></thead><tbody><tr><td>1/30</td><td>1166</td><td>1403</td><td>1.20</td></tr><tr><td>1/40</td><td>2056</td><td>2576</td><td>1.25</td></tr><tr><td>1/50</td><td>3337</td><td>4085</td><td>1.22</td></tr><tr><td>1/60</td><td>4710</td><td>5824</td><td>1.23</td></tr><tr><td>1/70</td><td>6638</td><td>8223</td><td>1.23</td></tr></tbody></table>

## 5.3 Conclusion and future directions

The objective of the present work is to propose reliable numerical schemes for the simulation of the Oldroyd-B model with Navier's slip boundary conditions. These formulations are based on a suitable factorization of the conformation tensor, which allows the derivation of an energy decay for the Oldroyd-B model. Such an estimate indicates, among other things, that the numerical schemes are stable. One of the main advantages of the schemes presented in this work is the possibility to adapt the time mesh with respect to $Re, We$, thus increasing the stability region. A stability analysis reveals that the size of the time step required for both schemes to decay exponentially is the same.

For the simulations presented in this article, two numerical experiments are performed: the first one numerically illustrates the decay of the energy $E_m$ with respect to time $t_m$, and the second one shows the effectiveness of the proposed schemes in a benchmark problem, namely the driven cavity.

The simulations demonstrate that both methods confirm the predictions of the theory. However, in terms of CPU time, the semi-implicit method appears to be superior to the method based on the Marchuk-Yanenko splitting algorithm. This is unusual, as in general, the splitting method is used to reduce CPU time. The main constraint in this task is to have a scheme that mimics the decay of the energy. Hence, the subdivisions in the splitting method presented are chosen accordingly. The current methodology should be further extended, particularly by: (i) devising an adapted discretization strategy, (ii) developing a competitive splitting approach that mimics the decay of the energy, (iii) devising a scalable solver, (iv) considering log-decomposition, which is possible since the conformation tensor is symmetric and positive definite, (v) developing an adapted discretization strategy for non-Newtonian fluids at high Weissemberg numbers, (vi) considering more realistic fluids, such as blood flow through a stenotic channel.

**conflicts of interests:** None

**References:**

[1] J.G. Abuga & T. Chinyoka. Benchmark solutions of the stabilized computations of flows of fluids governed by the Rolie-Poly constitutive model, *J. Phys. Commun.*, **4**, No 1, 2020.

[2] F. P.T. Baaijens. Mixed finite element methods for viscoelastic flow analysis: a review, *J. Non-Newtonian Fluid Mech.*, **79**, 361–385, 1998.

[3] J. Baranger, C. Guillope, & J. C. Saut. Mathematical analysis of differential models for viscoelastic fluids, in: J. M. Piau, J. F. Agassant (Eds). *Rheology for polymer melt processing*, Elsevier, Amsterdam, 199–236, 1996.

[4] J. Baranger, & D. Sandri. finite element approximation of viscoelastic fluid flow, *Numer. Math.*, **63**, 13–27, 1992.

[5] J. Barrett, & E. Suli. Existence of global weak solutions to dumbbell models for dilute polymers with microscopic cut-off, *M3AS.*, **18**, 935–971, 2008.

[6] M. Behr, L. Franca, & T. Tezduyar. Stabilized finite element methods for the velocity-pressure-stress formulation of incompressible flows, *Comput. Methods Appl. Mech. Engrg.*, **104**, 31–48, 1993.
[7] D. Boffi, F. Brezzi, & M. Fortin. Mixed Finite Element Methods and Applications, Springer Series in Computational Mathematics, Springer Verlag Berlin, 2013.

[8] S. Boyaval, T. Lelievre, & C. Mangoubi. Free energy dissipative schemes for the Oldroyd-B model, *ESAIM: M2AN*, DOI:10.1051/m2an/2009008, 2009.

[9] J.Y. Chemin, & N. Masmoudi. About lifespan of regular solutions of equations related to viscoelastic fluids, *SIAM J. Math. Anal.*, **33**, 84–112, 2001.

[10] P. G. Ciarlet. Three dimensional elasticity, *North Holland*, 1988.

[11] P. Ciarlet. Finite Element Method for Elliptic Problems, *North Holland, Amstedam*, 1978.

[12] C.A. Desoer & M. Vidyasagar. Feedback Systems Input-Output Properties, *Electrical Sciences. Academic Press, New York*, 1975.

[13] J.K. Djoko, J. Koko, & T. Sayah. Energy estimate for Oldroyd-B model under Tresca boundary condition: scheme preserving properties, *Comp. Appl. Math.*, **43**, 422, (2024).

[14] R. Fattal, & R. Kupferman. constitutive laws for the matrix-logarithm of the conformation tensor, *J. Non-Newtonian Fluid Mech.*, **123**, 281–285, 2005.

[15] M. Fortin, & D. Esselaoui. A finite element procedure for viscoelastic flows, *Int. J. Numer. Methods Fluids.*, **7**, 1035–1052, 1987.

[16] V. Girault, & P.A. Raviart. Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms, Springer-Verlag, Berlin, Heidelberg, New-York, Tokyo, 1986.

[17] A. Groisman, & V. Steinberg. Mechanism of elastic instability in Couette flow of polymer solutions, *Experiment Phys. Fluids.*, **10**, 2451–2463, 1998.

[18] M. A. Hulsen, M. A. Spanjaards, P. D. Anderson. Numerical stability of four positive(semi-)definite reformulations for viscoelastic fluid models in benchmark flows, *J. Non-Newtonian Fluid Mech.*, **297**, 104666, 2021.

[19] A. Lozinski, & R. G. Owens. An energy estimate for the Oldroyd B model: theory and applications, *J. Non-Newtonian Fluid Mech.*, **112**, 161–176, 2003.

[20] J. M. Marchal, & M. J. Crochet. A new mixed finite element for calculating viscoelastic flows, *J. Non-Newtonian Fluid Mech.*, **26**, 77–114, 1987.

[21] R. Fattal, & R. Kupferman. Time dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation, *J. Non-Newtonian Fluid Mech.*, **126**, 23–37, 2005.

[22] M. Fortin, & A. Fortin. A new approach for the FEM simulations of viscoelastic flows, *J. Non-Newtonian Fluid Mech.*, **32**, 295–310, 1989.

[23] R. Glowinski. Finite element methods for incompressible viscous flow, *In Handbook of Numerical Analysis, P.G Ciarlet and J.L Lions, eds, IX*, North Holland; Amsterdam, 3–1176, 2003.

[24] R. Glowinski, T.-W. Pan & X.-C. Tai. Some facts about operatorsplitting and alternating direction methods, in *Splitting Methods in Communication, Imaging, Science, and Engineering, Springer*, 19–94, 2016.

[25] F. Hecht. New development in FreeFem++, *J. Numer. Math.*, **20**, 251-266, 2012.

[26] P. Kumar. Displacement of Fluids in Permeable Media. *Wseas Transaction on Heat and Mass Transfer*, **17**, 2022.

[27] N. Kural, S. Shokhista, & A. Nabijon. Unsteady flow of rheologically complex fluids in cylindrical pipes, *Earth Sciences and Human Constructions*, **4**, 2024.

[28] P.L. Lions, & N. Masmoudi. Global solutions for some Oldroyd models of non-Newtonian flows, *Chinese Ann. Math.*, Ser. B. **21**, 131–146, 2000.

[29] F.H. Lin, C. Liu, & P. Zhang. On a micro-macro model for polymeric fluids near equilibrium, *Commun. Pure. Appl. Math.*, LX, 838–866, 2007.

[30] R. G. Owens. Computational Rheology. *Imperial College Press, London*, 2002.

[31] T.W. Pan, J. Hao & R. Glowinski. On the simulation of a time dependent cavity flow of an Oldroyd-B fluid, *Int. J. Numer. Fluids.*, **60**, 791–808, 2009.

[32] T.W. Pan, J. Hao & R. Glowinski. Positive Definiteness Preserving Approaches for Viscoelastic Flow of Oldroyd-B Fluids:
Applications to a Lid-Driven Cavity Flow and
a Particulate Flow, *Numerical Methods for
Non-Newtonian Fluids, Special Volume (R.
Glowinski and Jinchao Xu, Guest Editors)
of Handbook of Numerical Analysis, XVI*,
433–481, 2011.

[33] K. R. Rajagopal. On the implicit constitutive theories. *J. Fluid. Mech.*, **550**, 243–249 2006.

[34] B. Thomases, & M. Shelley. A transition to mixing and oscillations in a Stokesian viscoelastic flow, *Phys Rev. Lett.*, **103**, 094501, 2009.

[35] C. Le Roux. On flows of viscoelastic fluids of Oldroyd type with wall slip, *J. Math. Fluid. Mech.*, **465**, 723–730, 2018.

[36] M. Vidyasagar. Nonlinear Systems Analysis, Prentice Hall, Englewood Cliffs New Jersey, 2nd edition, 1993.

[37] R. Verfurth. A Review of A Posteriori Error Estimation and Adaptive Mesh-Refinement Techniques, *Wiley and Teubner Mathematics*, (1996).

[38] Y.J. Lee, J. Xu C. S. Zhang. Stable Finite Element Discretizations for Viscoelastic Flow Models, *In Handbook of Numerical Analysis, Numerical Methods for Non-Newtonian Fluids, eds, R. Glowinski, J. Xu., XVI*, North Holland; Amsterdam, 371–432, 2011.

[39] A. Srivastava, M. Bharty, & H. Mahato. Investigation of magneto-convection in viscoelastic fluid saturated anisotropic porous layer under local thermal non-equilibrium condition, *Int. J. of Applied Mathematics, Computational Science and Systems Engineering*, **6**, 2024.

**Contribution of Individual Authors to the**
**Creation of a Scientific Article (Ghostwriting**
**Policy)**

The authors equally contributed in the present research, at all stages from the formulation of the problem to the final findings and solution.

**Sources of Funding for Research Presented in a Scientific Article or Scientific Article Itself**

No funding was received for conducting this study.

**Conflict of Interest**

The authors have no conflicts of interest to declare that are relevant to the content of this article.

**Creative Commons Attribution License 4.0**
**(Attribution 4.0 International, CC BY 4.0)**

This article is published under the terms of the Creative Commons Attribution License 4.0
https://creativecommons.org/licenses/by/4.0/deed.en
_US