<figure><img src="image_26.png" /><figcaption>FIGURE 20. Velocity magnitude, pressure and $\operatorname{tr}(\boldsymbol{c})$ (in log-scale) fields. Flow past a confined cylinder, Oldroyd-B model, $Wi = 0.7$. Top to bottom: velocity magnitude, pressure and trace of the conformation tensor (in log-scale).</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>FIGURE 21. Profile of first normal polymer stress $\tau_{xx}$ along the curve coordinate $s$. Flow past a confined cylinder, Oldroyd-B model, $Wi = 0.6$. Left: profile along the cylinder surface and in the wake of the cylinder. Right: the same profile zoomed in the wake of the cylinder. (a) Profile along the curve coordinate. (b) Profile zoomed in the wake of the cylinder.</figcaption></figure>

### 5.2.1. Oldroyd-B model

We present the results for two Weissenberg numbers, $Wi = 0.6$ and $Wi = 0.7$. In both cases, the overall characteristics of the flow are the same so, for two-dimensional plots, we restrict ourselves to the two-dimensional fields (velocity, pressure and trace of the conformation tensor) on the MC4 mesh for the highest Weissenberg number, once the flow has reached the stationary state (Fig. 20). For $Wi = 0.7$, this occurs at time $t \approx 3$. We observe three regions of large polymer stresses. As discussed previously, a filament is present in the wake of the cylinder. Two large elongation zones also appear between the upper part of the cylinder and the opposite solid surface, at the neighborhood of the boundary, on both sides.
<figure><img src="image_28.png" /><figcaption>FIGURE 22. First normal polymer stress $\tau_{xx}$ profile zoomed in the wake of the cylinder. Flow past a confined cylinder, Oldroyd-B model, $Wi = 0.7$. Initial meshes (**left**) and refined meshes (**right**). (a) MC-type meshes. (b) MR-type meshes.</figcaption></figure>

To go further, we compare our results with those proposed in [20,32]. For $Wi = 0.6$, the stationary regime is obtained at time $t \approx 2$. Figure 21 shows the $xx$-component of the polymer stress,

$$ \tau_{xx} = \frac{1 - \beta}{Wi}(c_{xx} - 1), $$

as a function of the curvilinear coordinate $s$ along the line composed of the cylinder boundary and the part of the symmetry axis located in the wake of the cylinder. The two stagnation points at the cylinder boundary correspond to $s = 0$ for the front one and $s = \pi$ for the rear one. As expected from the two-dimensional plot, we observe two peaks for $\tau_{xx}$. The first one, more pronounced, is approximately centered at the throat between the cylinder and the solid surface, $s \approx \pi/2$, and the second one is located in the wake close to the cylinder. These results are in very good agreement with both [20,32].

We now turn to the case $Wi = 0.7$, for which a lack of convergence on $\tau_{xx}$ in the wake of the cylinder is generally reported in the literature (see [20,32]). This is observed in Figure 22 (**left**) with the MC-type meshes. We find a profile close to the one proposed in [32] with the MC1 mesh and close to the one proposed in [20] with the MC3 one. Beyond MC3, we need finer meshes that make it possible to capture thinner filaments with larger values of the conformation tensor. To explore this, we consider the second family of refined meshes, MR1, ..., MR7. The stress profiles are shown in Figure 22 (**right**), making it obvious that mesh convergence is very difficult to obtain and that the peak value increases as the mesh is refined. However, convergence seems close with the MR7 mesh.

### 5.2.2. FENE-CR model

Finally, we consider the FENE-CR fluid at $Wi = 0.7$ for different values of the elongation parameter $b$. Figure 23 shows the trace of the conformation tensor in log-scale. As for the lid-driven cavity, when the value of $b$ increases, the zones of high elongation become thinner and profiles look qualitatively similar to those obtained with the Oldroyd-B model.

Results obtained for the $xx$-component of the polymer stress,

$$ \tau_{xx} = \frac{1 - \beta}{Wi} \frac{b}{b - \text{tr}(c)} (c_{xx} - 1), $$
<figure><img src="image_28.png" /><figcaption>FIGURE 22. First normal polymer stress $\tau_{xx}$ profile zoomed in the wake of the cylinder. Flow past a confined cylinder, Oldroyd-B model, $Wi = 0.7$. Initial meshes (**left**) and refined meshes (**right**). (a) MC-type meshes. (b) MR-type meshes.</figcaption></figure>

To go further, we compare our results with those proposed in [20,32]. For $Wi = 0.6$, the stationary regime is obtained at time $t \approx 2$. Figure 21 shows the $xx$-component of the polymer stress,

$$ \tau_{xx} = \frac{1 - \beta}{Wi}(c_{xx} - 1), $$

as a function of the curvilinear coordinate $s$ along the line composed of the cylinder boundary and the part of the symmetry axis located in the wake of the cylinder. The two stagnation points at the cylinder boundary correspond to $s = 0$ for the front one and $s = \pi$ for the rear one. As expected from the two-dimensional plot, we observe two peaks for $\tau_{xx}$. The first one, more pronounced, is approximately centered at the throat between the cylinder and the solid surface, $s \approx \pi/2$, and the second one is located in the wake close to the cylinder. These results are in very good agreement with both [20,32].

We now turn to the case $Wi = 0.7$, for which a lack of convergence on $\tau_{xx}$ in the wake of the cylinder is generally reported in the literature (see [20,32]). This is observed in Figure 22 (**left**) with the MC-type meshes. We find a profile close to the one proposed in [32] with the MC1 mesh and close to the one proposed in [20] with the MC3 one. Beyond MC3, we need finer meshes that make it possible to capture thinner filaments with larger values of the conformation tensor. To explore this, we consider the second family of refined meshes, MR1, ..., MR7. The stress profiles are shown in Figure 22 (**right**), making it obvious that mesh convergence is very difficult to obtain and that the peak value increases as the mesh is refined. However, convergence seems close with the MR7 mesh.

### 5.2.2. FENE-CR model

Finally, we consider the FENE-CR fluid at $Wi = 0.7$ for different values of the elongation parameter $b$. Figure 23 shows the trace of the conformation tensor in log-scale. As for the lid-driven cavity, when the value of $b$ increases, the zones of high elongation become thinner and profiles look qualitatively similar to those obtained with the Oldroyd-B model.

Results obtained for the $xx$-component of the polymer stress,

$$ \tau_{xx} = \frac{1 - \beta}{Wi} \frac{b}{b - \text{tr}(c)} (c_{xx} - 1), $$
<figure><img src="image_29.png" /><figcaption>FIGURE 23. Trace of the conformation tensor (in log-scale). Flow past a confined cylinder, FENE-CR model, $Wi = 0.7$. Top to bottom: $b = 10$, $b = 100$ and $b = 1000$.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>FIGURE 24. First normal polymer stress $\tau_{xx}$ profile zoomed in near wake of the cylinder. Flow past a confined cylinder, FENE-CR model, $Wi = 0.7$. (a) $b = 100$. (b) $b = 1000$. Results obtained with the MC-type meshes.</figcaption></figure>

in the wake of the cylinder are shown in Figure 24, for $b = 100$ and $b = 1000$. For $b = 100$, mesh convergence is reached relatively easily, whereas convergence becomes slower for $b = 1000$. In this last case, difficulties are similar to those encountered with the Oldroyd-B model, with in fact similar values of the maximum value of the trace (about 187 for the Oldroyd-B model and 135 for the FENE-CR model). This is explained by the fact that the limit $b = 1000$ is far beyond the maximum value taken by $\text{tr}(c)$ with the Oldroyd-B model, so the prefactor $b/(b - \text{tr}(c))$ in the FENE-CR constitutive law remains close to one everywhere in the flow (contrary to what happens in the lid-driven cavity case).

## APPENDIX A. POSITIVITY OF A (DISCRETELY) TRANSPORTED TENSOR

The aim of this section is to prove that, with a first-order backward upwind scheme for the transport term, the solution of the transport equation for the conformation tensor (21d) remains definite positive, under a constraint for the time-step. We proceed in two steps. We first deal with the homogeneous transport (i.e. the
<figure><img src="image_29.png" /><figcaption>FIGURE 23. Trace of the conformation tensor (in log-scale). Flow past a confined cylinder, FENE-CR model, $Wi = 0.7$. Top to bottom: $b = 10$, $b = 100$ and $b = 1000$.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>FIGURE 24. First normal polymer stress $\tau_{xx}$ profile zoomed in near wake of the cylinder. Flow past a confined cylinder, FENE-CR model, $Wi = 0.7$. (a) $b = 100$. (b) $b = 1000$. Results obtained with the MC-type meshes.</figcaption></figure>

in the wake of the cylinder are shown in Figure 24, for $b = 100$ and $b = 1000$. For $b = 100$, mesh convergence is reached relatively easily, whereas convergence becomes slower for $b = 1000$. In this last case, difficulties are similar to those encountered with the Oldroyd-B model, with in fact similar values of the maximum value of the trace (about 187 for the Oldroyd-B model and 135 for the FENE-CR model). This is explained by the fact that the limit $b = 1000$ is far beyond the maximum value taken by $\text{tr}(c)$ with the Oldroyd-B model, so the prefactor $b/(b - \text{tr}(c))$ in the FENE-CR constitutive law remains close to one everywhere in the flow (contrary to what happens in the lid-driven cavity case).

## APPENDIX A. POSITIVITY OF A (DISCRETELY) TRANSPORTED TENSOR

The aim of this section is to prove that, with a first-order backward upwind scheme for the transport term, the solution of the transport equation for the conformation tensor (21d) remains definite positive, under a constraint for the time-step. We proceed in two steps. We first deal with the homogeneous transport (i.e. the
sub-problem obtained by setting the reaction terms to zero), then with the pure reaction equation (i.e. the sub-problem obtained by setting the convection velocity to zero), and finally merge the arguments to obtain the positivity result for the original equation.

We suppose that the fraction $g(c)/\lambda$ is equal to a constant, which we denote by $1/\lambda$ for short. However, the following arguments only need that $g(c)/\lambda$ be a known and non-negative value in each cell of the mesh. Hence, they also hold for the FENE-CR model if the time-discretization is such that the coefficient $g(c)/\lambda$ is computed with the beginning-of-step value of $c$, $c^n$, provided that $\text{tr}(c^n)$ is lower than $b$ in any cell of the mesh. This is indeed the case with the scheme used in this paper, at the price of constraints for the time step (see Sect. B), which will need to be satisfied as the same time as those appearing in this section.

## A.1. Pure transport

We first recall the arguments showing the positivity of the conformation tensor in the continuous setting, before mimicking the proof in the discrete setting.

Let $c(x, t) \in \mathbb{R}^{d \times d}$ be a tensor satisfying

$$ \partial_t c + \mathbf{u} \cdot \nabla c = 0 \quad \text{i.e.} \quad \partial_t c_{i,j} + \mathbf{u} \cdot \nabla c_{i,j} = 0 \text{ for } 1 \le i, j \le d, \qquad (A.1) $$

with $c(x, 0)$ positive and definite. The problem is supposed to be posed over a finite domain $\Omega \subset \mathbb{R}^d$ and, for the sake of simplicity, we suppose that $\mathbf{u} \cdot \mathbf{n}_{\partial\Omega} = 0$ over the whole domain boundary and that the velocity field $\mathbf{u}$ is divergence-free. Let $\mathbf{X} \in \mathbb{R}^d$ be a (fixed) vector. Then, for $1 \le i, j \le d$, we have

$$ X_i X_j [\partial_t c_{i,j} + \mathbf{u} \cdot \nabla c_{i,j}] = \partial_t (c_{i,j} X_i X_j) + \mathbf{u} \cdot \nabla (c_{i,j} X_i X_j) = 0. $$

Summing over the indices $i$ and $j$ yields

$$ \partial_t((c\mathbf{X}, \mathbf{X})) + \mathbf{u} \cdot \nabla((c\mathbf{X}, \mathbf{X})) = 0. $$

Since the transport equation satisfies a maximum principle and that $(c(x, 0)\mathbf{X}, \mathbf{X}) > 0$ over $\Omega$ by assumption on the initial data, we get that $(c(x, t)\mathbf{X}, \mathbf{X}) > 0$ over $\Omega \times \mathbb{R}$, i.e. that the tensor $c$ remains positive definite at all locations and times. The following lemma states a discrete counterpart to this property.

**Lemma A.1.** The discretization of (A.1) by the first-order backward-Euler upwind scheme preserves the positivity of the tensor $c$.

*Proof.* The discretization of (A.1) by the first-order backward-Euler upwind scheme reads

$$ \text{for } 1 \le i, j \le d, \quad \forall K \in \mathcal{M}, \quad \frac{|K|}{\delta t} \left( (c_{i,j})_K^{n+1} - (c_{i,j})_K^n \right) + \sum_{\sigma=K|L} |\sigma| (c_{i,j})_{\sigma}^{n+1} u_{K,\sigma} = 0, $$

where $u_{K,\sigma}$ is an approximation of $\mathbf{u} \cdot \mathbf{n}_{K,\sigma}$ such that the following divergence-free constraint is satisfied (remember that $\mathbf{u}$ is supposed to be divergence-free):

$$ \sum_{\sigma=K|L} |\sigma| u_{K,\sigma} = 0. \qquad (A.2) $$

and

$$ (c_{i,j})_{\sigma}^{n+1} = (c_{i,j})_{K}^{n+1} \quad \text{if } u_{K,\sigma} \ge 0, \quad (c_{i,j})_{\sigma}^{n+1} = (c_{i,j})_{L}^{n+1} \quad \text{otherwise.} \qquad (A.3) $$

Note that, to simplify notations, we do not associate a time-index to the velocity, but the fact that the velocity depends on time is indifferent in the present proof. Let $\mathbf{X} \in \mathbb{R}^d$ be a (fixed) vector. Then, for $K \in \mathcal{M}$ and $1 \le i, j \le d$, we have

$$ \frac{|K|}{\delta t} \left( X_i X_j (c_{i,j})_K^{n+1} - X_i X_j (c_{i,j})_K^n \right) + \sum_{\sigma=K|L} |\sigma| X_i X_j (c_{i,j})_{\sigma}^{n+1} u_{K,\sigma} = 0. $$
sub-problem obtained by setting the reaction terms to zero), then with the pure reaction equation (i.e. the sub-problem obtained by setting the convection velocity to zero), and finally merge the arguments to obtain the positivity result for the original equation.

We suppose that the fraction $g(c)/\lambda$ is equal to a constant, which we denote by $1/\lambda$ for short. However, the following arguments only need that $g(c)/\lambda$ be a known and non-negative value in each cell of the mesh. Hence, they also hold for the FENE-CR model if the time-discretization is such that the coefficient $g(c)/\lambda$ is computed with the beginning-of-step value of $c$, $c^n$, provided that tr($c^n$) is lower than $b$ in any cell of the mesh. This is indeed the case with the scheme used in this paper, at the price of constraints for the time step (see Sect. B), which will need to be satisfied as the same time as those appearing in this section.

## A.1. Pure transport

We first recall the arguments showing the positivity of the conformation tensor in the continuous setting, before mimicking the proof in the discrete setting.

Let $c(x, t) \in \mathbb{R}^{d \times d}$ be a tensor satisfying

$$ \partial_t c + \mathbf{u} \cdot \nabla c = 0 \quad \text{i.e.} \quad \partial_t c_{i,j} + \mathbf{u} \cdot \nabla c_{i,j} = 0 \text{ for } 1 \le i, j \le d, \qquad (A.1) $$

with $c(x, 0)$ positive and definite. The problem is supposed to be posed over a finite domain $\Omega \subset \mathbb{R}^d$ and, for the sake of simplicity, we suppose that $\mathbf{u} \cdot \mathbf{n}_{\partial\Omega} = 0$ over the whole domain boundary and that the velocity field $\mathbf{u}$ is divergence-free. Let $\mathbf{X} \in \mathbb{R}^d$ be a (fixed) vector. Then, for $1 \le i, j \le d$, we have

$$ X_i X_j [\partial_t c_{i,j} + \mathbf{u} \cdot \nabla c_{i,j}] = \partial_t (c_{i,j} X_i X_j) + \mathbf{u} \cdot \nabla (c_{i,j} X_i X_j) = 0. $$

Summing over the indices $i$ and $j$ yields

$$ \partial_t ((c\mathbf{X}, \mathbf{X})) + \mathbf{u} \cdot \nabla ((c\mathbf{X}, \mathbf{X})) = 0. $$

Since the transport equation satisfies a maximum principle and that $(c(x, 0)\mathbf{X}, \mathbf{X}) > 0$ over $\Omega$ by assumption on the initial data, we get that $(c(x, t)\mathbf{X}, \mathbf{X}) > 0$ over $\Omega \times \mathbb{R}$, i.e. that the tensor $c$ remains positive definite at all locations and times. The following lemma states a discrete counterpart to this property.

**Lemma A.1.** The discretization of (A.1) by the first-order backward-Euler upwind scheme preserves the positivity of the tensor $c$.

*Proof.* The discretization of (A.1) by the first-order backward-Euler upwind scheme reads

$$ \text{for } 1 \le i, j \le d, \quad \forall K \in \mathcal{M}, \quad \frac{|K|}{\delta t} \left( (c_{i,j})_K^{n+1} - (c_{i,j})_K^n \right) + \sum_{\sigma=K|L} |\sigma| (c_{i,j})_{\sigma}^{n+1} u_{K,\sigma} = 0, $$

where $u_{K,\sigma}$ is an approximation of $\mathbf{u} \cdot \mathbf{n}_{K,\sigma}$ such that the following divergence-free constraint is satisfied (remember that $\mathbf{u}$ is supposed to be divergence-free):

$$ \sum_{\sigma=K|L} |\sigma| u_{K,\sigma} = 0. \qquad (A.2) $$

and

$$ (c_{i,j})_{\sigma}^{n+1} = (c_{i,j})_{K}^{n+1} \text{ if } u_{K,\sigma} \ge 0, \quad (c_{i,j})_{\sigma}^{n+1} = (c_{i,j})_{L}^{n+1} \text{ otherwise.} \qquad (A.3) $$

Note that, to simplify notations, we do not associate a time-index to the velocity, but the fact that the velocity depends on time is indifferent in the present proof. Let $\mathbf{X} \in \mathbb{R}^d$ be a (fixed) vector. Then, for $K \in \mathcal{M}$ and $1 \le i, j \le d$, we have

$$ \frac{|K|}{\delta t} \left( X_i X_j (c_{i,j})_K^{n+1} - X_i X_j (c_{i,j})_K^n \right) + \sum_{\sigma=K|L} |\sigma| X_i X_j (c_{i,j})_{\sigma}^{n+1} u_{K,\sigma} = 0. $$
Let us define $S_K^n = (c_K^n X, X)$. Then, summing the previous relation over $i$ and $j$ yields

$$ \frac{|K|}{\delta t} (S_K^{n+1} - S_K^n) + \sum_{\sigma=K|L} |\sigma| S_{\sigma}^{n+1} u_{K,\sigma} = 0, $$

with, thanks to definition (A.3), $S_{\sigma}^{n+1} = S_K^{n+1}$ if $u_{K,\sigma} \ge 0$ and $S_{\sigma}^{n+1} = S_L^{n+1}$ otherwise. We may see $S$ as an unknown obeying a first-order backward-Euler upwind scheme, which is known to preserve the positivity. For any reasonable initialization, we have $S_K^0 > 0$ for $K \in \mathcal{M}$, so $S_K^n > 0$ for $n > 0$, which means that the discrete tensor $c$ is positive in each cell and at any time. $\square$

**Remark A.2.** In the case with the change of variable to the log-conformation tensor, taking the exponential after the advection step allows to automatically ensure the positivity of $c$.

## A.2. The pure reaction problem

With the implicit variant, the ODE step for the conformation tensor can be written in each cell, dropping for short the cell index,

$$ \frac{1}{\delta t} (\mathbf{c}^{n+1} - \mathbf{c}^n) = \nabla \mathbf{u} \mathbf{c}^{n+1} + (\mathbf{c}^{n+1})^t (\nabla \mathbf{u})^t - \frac{1}{\lambda} (\mathbf{c}^{n+1} - \mathbf{I}_d). \quad (A.4) $$

Note that, in the second term on the right-hand side, we have replaced $\mathbf{c}^{n+1}$ by $(\mathbf{c}^{n+1})^t$. This is justified by the fact that, provided that $\mathbf{c}^n$ is symmetrical, any solution of (A.4) is also symmetrical and thus satisfies the initial equation. Equation (A.4) is a linear system for the components of $\mathbf{c}^{n+1}$ which is regular for $\delta t$ small enough. For instance, the condition

$$ \frac{1}{\delta t} + \frac{1}{\lambda} > 2 \max_{1 \le i \le d} \sum_{j=1}^{d} \left| ((\nabla \mathbf{u})_{i,j}) \right| \quad (A.5) $$

ensures that the diagonal of the matrix of this system is strictly dominant. The following lemma states that the scheme (A.4) preserves the definite positivity of the unknown.

**Lemma A.3.** Let us suppose that the time step is such that Condition (A.5) is satisfied, and that $\mathbf{c}^n$ is a symmetric positive definite tensor. Then so is $\mathbf{c}^{n+1}$.

*Proof.* Let $\alpha \mapsto \mathbf{c}(\alpha)$, $[0, 1] \to \mathbb{R}^{d \times d}$ be the function defined by taking $\mathbf{c}(\alpha)$ as the solution to equation (A.4) with $\alpha \delta t$ instead of $\delta t$:

$$ \frac{1}{\delta t} (\mathbf{c}(\alpha) - \mathbf{c}^n) = \alpha \left[ \nabla \mathbf{u} \mathbf{c}(\alpha) + \mathbf{c}(\alpha)^t (\nabla \mathbf{u})^t - \frac{1}{\lambda} (\mathbf{c}(\alpha) - \mathbf{I}_d) \right]. \quad (A.6) $$

Since, thanks to condition (A.5), $\mathbf{c}(\alpha)$ is the solution to a regular linear system, the function $\mathbf{c}(\alpha)$ is regular. In addition, $\mathbf{c}(0) = \mathbf{c}^n$ is a symmetric positive definite tensor. Let us suppose that $\mathbf{c}(\alpha)$ is not definite positive for some $\alpha \in [0, 1]$, which means that one at least of its eigenvalues is non-positive. Then, by continuity of $\alpha \mapsto \min(\text{eigenvalues of } \mathbf{c}(\alpha))$, there exists $\alpha_0 \in (0, 1]$ such that $\mathbf{c}(\alpha_0)$ has a zero eigenvalue. Let $\mathbf{X} \in \mathbb{R}^d$ be a unit eigenvector associated with this zero eigenvalue. We have, by equation (A.6),

$$ \begin{aligned} -\frac{1}{\delta t} (\mathbf{c}^n \mathbf{X}, \mathbf{X}) &= \alpha_0 \left[ (\nabla \mathbf{u} \mathbf{c}(\alpha_0) \mathbf{X}, \mathbf{X}) + (\mathbf{c}(\alpha_0)^t (\nabla \mathbf{u})^t \mathbf{X}, \mathbf{X}) + \frac{1}{\lambda} \right] \\ &= \alpha_0 \left( 2 (\nabla \mathbf{u} \mathbf{c}(\alpha_0) \mathbf{X}, \mathbf{X}) + \frac{1}{\lambda} \right) = \frac{\alpha_0}{\lambda}, \end{aligned} $$

which raises a contradiction since, by assumption, $(\mathbf{c}^n \mathbf{X}, \mathbf{X}) > 0$. $\square$
Let us define $S_K^n = (c_K^n X, X)$. Then, summing the previous relation over $i$ and $j$ yields

$$ \frac{|K|}{\delta t} (S_K^{n+1} - S_K^n) + \sum_{\sigma=K|L} |\sigma| S_{\sigma}^{n+1} u_{K,\sigma} = 0, $$

with, thanks to definition (A.3), $S_{\sigma}^{n+1} = S_K^{n+1}$ if $u_{K,\sigma} \ge 0$ and $S_{\sigma}^{n+1} = S_L^{n+1}$ otherwise. We may see $S$ as an unknown obeying a first-order backward-Euler upwind scheme, which is known to preserve the positivity. For any reasonable initialization, we have $S_K^0 > 0$ for $K \in \mathcal{M}$, so $S_K^n > 0$ for $n > 0$, which means that the discrete tensor $c$ is positive in each cell and at any time. □

**Remark A.2.** In the case with the change of variable to the log-conformation tensor, taking the exponential after the advection step allows to automatically ensure the positivity of $c$.

## A.2. The pure reaction problem

With the implicit variant, the ODE step for the conformation tensor can be written in each cell, dropping for short the cell index,

$$ \frac{1}{\delta t} (\mathbf{c}^{n+1} - \mathbf{c}^n) = \nabla \mathbf{u} \mathbf{c}^{n+1} + (\mathbf{c}^{n+1})^t (\nabla \mathbf{u})^t - \frac{1}{\lambda} (\mathbf{c}^{n+1} - \mathbf{I}_d). \quad (A.4) $$

Note that, in the second term on the right-hand side, we have replaced $\mathbf{c}^{n+1}$ by $(\mathbf{c}^{n+1})^t$. This is justified by the fact that, provided that $\mathbf{c}^n$ is symmetrical, any solution of (A.4) is also symmetrical and thus satisfies the initial equation. Equation (A.4) is a linear system for the components of $\mathbf{c}^{n+1}$ which is regular for $\delta t$ small enough. For instance, the condition

$$ \frac{1}{\delta t} + \frac{1}{\lambda} > 2 \max_{1 \le i \le d} \sum_{j=1}^{d} \left| \left( (\nabla \mathbf{u})_{i,j} \right) \right| \quad (A.5) $$

ensures that the diagonal of the matrix of this system is strictly dominant. The following lemma states that the scheme (A.4) preserves the definite positivity of the unknown.

**Lemma A.3.** Let us suppose that the time step is such that Condition (A.5) is satisfied, and that $\mathbf{c}^n$ is a symmetric positive definite tensor. Then so is $\mathbf{c}^{n+1}$.

*Proof.* Let $\alpha \mapsto \mathbf{c}(\alpha)$, $[0, 1] \to \mathbb{R}^{d \times d}$ be the function defined by taking $\mathbf{c}(\alpha)$ as the solution to equation (A.4) with $\alpha \delta t$ instead of $\delta t$:

$$ \frac{1}{\delta t} (\mathbf{c}(\alpha) - \mathbf{c}^n) = \alpha \left[ \nabla \mathbf{u} \mathbf{c}(\alpha) + \mathbf{c}(\alpha)^t (\nabla \mathbf{u})^t - \frac{1}{\lambda} (\mathbf{c}(\alpha) - \mathbf{I}_d) \right]. \quad (A.6) $$

Since, thanks to condition (A.5), $\mathbf{c}(\alpha)$ is the solution to a regular linear system, the function $\mathbf{c}(\alpha)$ is regular. In addition, $\mathbf{c}(0) = \mathbf{c}^n$ is a symmetric positive definite tensor. Let us suppose that $\mathbf{c}(\alpha)$ is not definite positive for some $\alpha \in [0, 1]$, which means that one at least of its eigenvalues is non-positive. Then, by continuity of $\alpha \mapsto \min(\text{eigenvalues of } \mathbf{c}(\alpha))$, there exists $\alpha_0 \in (0, 1]$ such that $\mathbf{c}(\alpha_0)$ has a zero eigenvalue. Let $\mathbf{X} \in \mathbb{R}^d$ be a unit eigenvector associated with this zero eigenvalue. We have, by equation (A.6),

$$ \begin{aligned} -\frac{1}{\delta t} (\mathbf{c}^n X, X) &= \alpha_0 \left[ (\nabla \mathbf{u} \mathbf{c}(\alpha_0) X, X) + (\mathbf{c}(\alpha_0)^t (\nabla \mathbf{u})^t X, X) + \frac{1}{\lambda} \right] \\ &= \alpha_0 \left( 2 (\nabla \mathbf{u} \mathbf{c}(\alpha_0) X, X) + \frac{1}{\lambda} \right) = \frac{\alpha_0}{\lambda}, \end{aligned} $$

which raises a contradiction since, by assumption, $(\mathbf{c}^n X, X) > 0$. □
### A.3. The complete transport-reaction equation

Let us now consider the complete transport equation for the conformation tensor:

$$ \forall K \in \mathcal{M}, \quad \frac{1}{\delta t} (c_K^{n+1} - c_K^n) + \frac{1}{|K|} \sum_{\sigma=K|L} |\sigma| c_{\sigma}^{n+1} u_{K,\sigma} = \nabla u_K c_K^{n+1} + (c_K^{n+1})^t (\nabla u_K)^t - \frac{1}{\lambda} (c_K^{n+1} - I_d), \quad (A.7) $$

where the velocity field is supposed to be discretely divergence-free (i.e. satisfies Eq. (A.2)) and the conformation tensor value at the face $\sigma$, $c_{\sigma}^{n+1}$, is given by the first-order backward upwind choice (A.3). We introduce the following set of constraints for the time step:

$$ \forall K \in \mathcal{M}, \quad \frac{1}{\delta t} + \frac{1}{\lambda} > 2 \max_{1 \le i \le d} \sum_{j=1}^{d} |(\nabla u_K)_{i,j}|. \quad (A.8) $$

We are now in position to state the following result.

**Lemma A.4.** Let the time step be such that equation (A.8) are verified, and let us assume that, for $K \in \mathcal{M}$, $c_K^n$ is a symmetric positive definite tensor. Then the linear System (A.7) is uniquely solvable and, for $K \in \mathcal{M}$, the tensor $c_K^{n+1}$ solution of (A.7) is symmetric positive definite.

*Proof.* The proof of this result is obtained by combining the arguments of the two previous sections. First, for each cell $K$, the upwind choice (A.3) implies that the convection operator produces, in the system for the tensor components associated with $K$, an additional non-negative term on the diagonal and non-positive off-diagonal additional entries. Moreover, thanks to relation (A.2), the diagonal additional term is equal to the sum of the off-diagonal ones. Inequality (A.8) thus still implies that the diagonal of the matrix of this system is dominant.

Then, as in Section A.2, let us multiply the right-hand side of (A.7) by $\alpha \in [0, 1]$; we easily check that the obtained system is still uniquely solvable, and denote by $(c_K^\alpha)_{K \in \mathcal{M}}$ its solution. By Section A.1, we know that, for $\alpha = 0$, $c_K^\alpha$ is symmetric definite positive, for any $K \in \mathcal{M}$. If, for $\alpha = 1$, $c_K^\alpha$ is not definite positive in at least one cell, by continuity, there exists $\alpha_0 \in (0, 1]$ and $K_0 \in \mathcal{M}$ such that $c_L^{\alpha_0}$ is non-negative in all cells $L \in \mathcal{M}$ and $c_{K_0}^{\alpha_0}$ has a zero eigenvalue. Let $X \in \mathbb{R}^d$ be a unit eigenvector associated with this zero eigenvalue. We have, by the same computations as in Sections A.2 and A.1,

$$ -\frac{1}{\delta t} (c_{K_0}^n X, X) + \frac{1}{|K|} \sum_{\sigma \in \mathcal{E}(K_0)} |\sigma| (c_{\sigma}^{\alpha_0} X, X) u_{K,\sigma} = \frac{\alpha_0}{\lambda}, $$

and, since $c_{K_0}^n$ is definite positive, $(c_{K_0}^{\alpha_0} X, X) = 0$ and $(c_L^{\alpha_0} X, X) \ge 0$ for $L \in \mathcal{M}$, we easily check that the upwind choice for $c_\sigma^{\alpha_0}$ implies that the left-hand side is negative, which raises a contradiction. $\square$

## APPENDIX B. STABILITY OF THE ODE STEP FOR THE FENE-CR MODEL

We study in this section the stability of the following semi-implicit scheme for the ODE step, for the FENE-CR model:

$$ \frac{1}{\delta t} (c^{n+1} - c^n) = \nabla u \ c^{n+1} + (c^{n+1})^t (\nabla u)^t - \frac{b}{\lambda(b - \text{tr}(c^n))} (c^{n+1} - I_d), \quad (B.9) $$

where $b$ is a real number strictly greater than the space dimension $d$. From Section A.2, we easily deduce that, provided that $c^n$ is a symmetric definite positive tensor such that $\text{tr}(c^n) < b$, then $c^{n+1}$ is also symmetric definite positive, under the stability constraint for the time step:

$$ \frac{1}{\delta t} > 2 \max_{1 \le i \le d} \sum_{j=1}^{d} |(\nabla u)_{i,j}|. \quad (B.10) $$

We search for an upper bound for $\delta t$ that ensures that $\text{tr}(c^{n+1}) < b$. To this purpose, we will use the following lemma.
### A.3. The complete transport-reaction equation

Let us now consider the complete transport equation for the conformation tensor:

$$ \forall K \in \mathcal{M}, \quad \frac{1}{\delta t} (c_K^{n+1} - c_K^n) + \frac{1}{|K|} \sum_{\sigma=K|L} |\sigma| c_{\sigma}^{n+1} u_{K,\sigma} = \nabla u_K c_K^{n+1} + (c_K^{n+1})^t (\nabla u_K)^t - \frac{1}{\lambda} (c_K^{n+1} - I_d), \quad (A.7) $$

where the velocity field is supposed to be discretely divergence-free (i.e. satisfies Eq. (A.2)) and the conformation tensor value at the face $\sigma$, $c_{\sigma}^{n+1}$, is given by the first-order backward upwind choice (A.3). We introduce the following set of constraints for the time step:

$$ \forall K \in \mathcal{M}, \quad \frac{1}{\delta t} + \frac{1}{\lambda} > 2 \max_{1 \le i \le d} \sum_{j=1}^{d} |(\nabla u_K)_{i,j}|. \quad (A.8) $$

We are now in position to state the following result.

**Lemma A.4.** Let the time step be such that equation (A.8) are verified, and let us assume that, for $K \in \mathcal{M}$, $c_K^n$ is a symmetric positive definite tensor. Then the linear System (A.7) is uniquely solvable and, for $K \in \mathcal{M}$, the tensor $c_K^{n+1}$ solution of (A.7) is symmetric positive definite.

*Proof.* The proof of this result is obtained by combining the arguments of the two previous sections. First, for each cell $K$, the upwind choice (A.3) implies that the convection operator produces, in the system for the tensor components associated with $K$, an additional non-negative term on the diagonal and non-positive off-diagonal additional entries. Moreover, thanks to relation (A.2), the diagonal additional term is equal to the sum of the off-diagonal ones. Inequality (A.8) thus still implies that the diagonal of the matrix of this system is dominant.

Then, as in Section A.2, let us multiply the right-hand side of (A.7) by $\alpha \in [0, 1]$; we easily check that the obtained system is still uniquely solvable, and denote by $(c_K^\alpha)_{K \in \mathcal{M}}$ its solution. By Section A.1, we know that, for $\alpha = 0$, $c_K^\alpha$ is symmetric definite positive, for any $K \in \mathcal{M}$. If, for $\alpha = 1$, $c_K^\alpha$ is not definite positive in at least one cell, by continuity, there exists $\alpha_0 \in (0, 1]$ and $K_0 \in \mathcal{M}$ such that $c_L^{\alpha_0}$ is non-negative in all cells $L \in \mathcal{M}$ and $c_{K_0}^{\alpha_0}$ has a zero eigenvalue. Let $X \in \mathbb{R}^d$ be a unit eigenvector associated with this zero eigenvalue. We have, by the same computations as in Sections A.2 and A.1,

$$ -\frac{1}{\delta t} (c_{K_0}^n X, X) + \frac{1}{|K|} \sum_{\sigma \in \mathcal{E}(K_0)} |\sigma| (c_{\sigma}^{\alpha_0} X, X) u_{K,\sigma} = \frac{\alpha_0}{\lambda}, $$

and, since $c_{K_0}^n$ is definite positive, $(c_{K_0}^{\alpha_0} X, X) = 0$ and $(c_L^{\alpha_0} X, X) \ge 0$ for $L \in \mathcal{M}$, we easily check that the upwind choice for $c_\sigma^{\alpha_0}$ implies that the left-hand side is negative, which raises a contradiction. $\square$

## APPENDIX B. STABILITY OF THE ODE STEP FOR THE FENE-CR MODEL

We study in this section the stability of the following semi-implicit scheme for the ODE step, for the FENE-CR model:

$$ \frac{1}{\delta t} (c^{n+1} - c^n) = \nabla u \ c^{n+1} + (c^{n+1})^t (\nabla u)^t - \frac{b}{\lambda(b - \text{tr}(c^n))} (c^{n+1} - I_d), \quad (B.9) $$

where $b$ is a real number strictly greater than the space dimension $d$. From Section A.2, we easily deduce that, provided that $c^n$ is a symmetric definite positive tensor such that $\text{tr}(c^n) < b$, then $c^{n+1}$ is also symmetric definite positive, under the stability constraint for the time step:

$$ \frac{1}{\delta t} > 2 \max_{1 \le i \le d} \sum_{j=1}^{d} |(\nabla u)_{i,j}|. \quad (B.10) $$

We search for an upper bound for $\delta t$ that ensures that $\text{tr}(c^{n+1}) < b$. To this purpose, we will use the following lemma.
**Lemma B.1.** Let $\mathbf{A}$ be a symmetric positive $d \times d$ matrix of entries $(A_{i,j})_{1\le i,j\le d}$. Then:

$$ ||\mathbf{A}||_{L^1} = \sum_{i,j=1}^{d} |A_{i,j}| \le d \operatorname{tr}(\mathbf{A}). $$

*Proof.* Let $\operatorname{sgn}(\mathbf{A})$ be the matrix of entries $(\operatorname{sgn}(A_{i,j}))_{1\le i,j\le d}$. Then $||\mathbf{A}||_{L^1} = \mathbf{A} : \operatorname{sgn}(\mathbf{A})$. In addition, $\operatorname{sgn}(\mathbf{A})$ is a symmetric diagonalizable matrix and its eigenvalues $\alpha_i$, $1 \le i \le d$, satisfy

$$ |\alpha_i| \le \max_{\mathbf{X} \in \mathbb{R}^d} \frac{|\operatorname{sgn}(\mathbf{A})\mathbf{X} \cdot \mathbf{X}|}{||\mathbf{X}||^2}. $$

We have, since the entries of $\operatorname{sgn}(\mathbf{A})$ are either 1 or -1,

$$ |\operatorname{sgn}(\mathbf{A})\mathbf{X} \cdot \mathbf{X}| = \left| \sum_{i=1}^{d} X_i \sum_{j=1}^{d} \operatorname{sgn}(A_{i,j}) X_j \right| \le \sum_{i,j=1}^{d} |X_i X_j| \le \frac{1}{2} \sum_{i,j=1}^{d} (X_i^2 + X_j^2) = d ||\mathbf{X}||^2. $$

Thus $|\alpha_i| \le d$, for $1 \le i \le d$. Since $\mathbf{A}$ is symmetric and positive, all its eigenvalues are non-negative and the results follows from the matrix Ky Fan inequality:

$$ \mathbf{A} : \operatorname{sgn}(\mathbf{A}) \le \Gamma(\mathbf{A}) \cdot \Gamma(\operatorname{sgn}(\mathbf{A})), $$

where, for a symmetric matrix $\mathbf{B}$, $\Gamma(\mathbf{B})$ stands for the vector obtained from the ordered set of the eigenvalues of $\mathbf{B}$. $\square$

Taking the trace of equation (B.9), we get

$$ \frac{1}{\delta t} (\operatorname{tr}(\mathbf{c}^{n+1}) - \operatorname{tr}(\mathbf{c}^n)) = 2 \nabla \mathbf{u} : \mathbf{c}^{n+1} - \frac{b}{\lambda(b - \operatorname{tr}(\mathbf{c}^n))} (\operatorname{tr}(\mathbf{c}^{n+1}) - d). $$

Under the constraint (B.10) for the time-step, the linear System (B.9) for $\mathbf{c}^{n+1}$ is regular. Thus, by continuity, since $\operatorname{tr}(\mathbf{c}^n) < b$, $\operatorname{tr}(\mathbf{c}^{n+1}) < b$ for $\delta t$ small enough. Let us suppose that there exists some values of the time step satisfying (B.10) and such that $\operatorname{tr}(\mathbf{c}^{n+1}) \ge b$; still by continuity, we may define $\delta t_0$ as the smallest time step such that $\operatorname{tr}(\mathbf{c}^{n+1}) = b$, which yields:

$$ \frac{1}{\delta t_0} (b - \operatorname{tr}(\mathbf{c}^n)) = 2 \nabla \mathbf{u} : \mathbf{c}^{n+1} - \frac{b}{\lambda(b - \operatorname{tr}(\mathbf{c}^n))} (b - d). $$

We thus have

$$ \frac{1}{\delta t_0} (b - \operatorname{tr}(\mathbf{c}^n)) + \frac{b(b-d)}{\lambda(b - \operatorname{tr}(\mathbf{c}^n))} \le 2 ||\nabla \mathbf{u}||_{L^\infty} ||\mathbf{c}^{n+1}||_{L^1}, $$

so, thanks to Lemma B.1,

$$ \frac{1}{\delta t_0} (b - \operatorname{tr}(\mathbf{c}^n)) + \frac{b(b-d)}{\lambda(b - \operatorname{tr}(\mathbf{c}^n))} \le 2d ||\nabla \mathbf{u}||_{L^\infty} \operatorname{tr}(\mathbf{c}^{n+1}) = 2db ||\nabla \mathbf{u}||_{L^\infty}. \quad (B.11) $$

We observe that the two terms at the left-hand side are positive. This inequality is thus impossible, and we have $\operatorname{tr}(\mathbf{c}^{n+1}) < b$ for any value of the time step, if

$$ \frac{b-d}{2d\lambda(b - \operatorname{tr}(\mathbf{c}^n))} > ||\nabla \mathbf{u}||_{L^\infty}. \quad (B.12) $$
**Lemma B.1.** Let $\mathbf{A}$ be a symmetric positive $d \times d$ matrix of entries $(A_{i,j})_{1\le i,j\le d}$. Then:

$$ \|\mathbf{A}\|_{L^1} = \sum_{i,j=1}^{d} |A_{i,j}| \le d \operatorname{tr}(\mathbf{A}). $$

*Proof.* Let $\operatorname{sgn}(\mathbf{A})$ be the matrix of entries $(\operatorname{sgn}(A_{i,j}))_{1\le i,j\le d}$. Then $\|\mathbf{A}\|_{L^1} = \mathbf{A} : \operatorname{sgn}(\mathbf{A})$. In addition, $\operatorname{sgn}(\mathbf{A})$ is a symmetric diagonalizable matrix and its eigenvalues $\alpha_i$, $1 \le i \le d$, satisfy

$$ |\alpha_i| \le \max_{\mathbf{X} \in \mathbb{R}^d} \frac{|\operatorname{sgn}(\mathbf{A})\mathbf{X} \cdot \mathbf{X}|}{\|\mathbf{X}\|^2}. $$

We have, since the entries of $\operatorname{sgn}(\mathbf{A})$ are either 1 or $-1$,

$$ |\operatorname{sgn}(\mathbf{A})\mathbf{X} \cdot \mathbf{X}| = \left| \sum_{i=1}^{d} X_i \sum_{j=1}^{d} \operatorname{sgn}(A_{i,j}) X_j \right| \le \sum_{i,j=1}^{d} |X_i X_j| \le \frac{1}{2} \sum_{i,j=1}^{d} (X_i^2 + X_j^2) = d \|\mathbf{X}\|^2. $$

Thus $|\alpha_i| \le d$, for $1 \le i \le d$. Since $\mathbf{A}$ is symmetric and positive, all its eigenvalues are non-negative and the results follows from the matrix Ky Fan inequality:

$$ \mathbf{A} : \operatorname{sgn}(\mathbf{A}) \le \Gamma(\mathbf{A}) \cdot \Gamma(\operatorname{sgn}(\mathbf{A})), $$

where, for a symmetric matrix $\mathbf{B}$, $\Gamma(\mathbf{B})$ stands for the vector obtained from the ordered set of the eigenvalues of $\mathbf{B}$. $\square$

Taking the trace of equation (B.9), we get

$$ \frac{1}{\delta t} (\operatorname{tr}(\mathbf{c}^{n+1}) - \operatorname{tr}(\mathbf{c}^n)) = 2 \nabla \mathbf{u} : \mathbf{c}^{n+1} - \frac{b}{\lambda(b - \operatorname{tr}(\mathbf{c}^n))} (\operatorname{tr}(\mathbf{c}^{n+1}) - d). $$

Under the constraint (B.10) for the time-step, the linear System (B.9) for $\mathbf{c}^{n+1}$ is regular. Thus, by continuity, since $\operatorname{tr}(\mathbf{c}^n) < b$, $\operatorname{tr}(\mathbf{c}^{n+1}) < b$ for $\delta t$ small enough. Let us suppose that there exists some values of the time step satisfying (B.10) and such that $\operatorname{tr}(\mathbf{c}^{n+1}) \ge b$; still by continuity, we may define $\delta t_0$ as the smallest time step such that $\operatorname{tr}(\mathbf{c}^{n+1}) = b$, which yields:

$$ \frac{1}{\delta t_0} (b - \operatorname{tr}(\mathbf{c}^n)) = 2 \nabla \mathbf{u} : \mathbf{c}^{n+1} - \frac{b}{\lambda(b - \operatorname{tr}(\mathbf{c}^n))} (b - d). $$

We thus have

$$ \frac{1}{\delta t_0} (b - \operatorname{tr}(\mathbf{c}^n)) + \frac{b(b-d)}{\lambda(b - \operatorname{tr}(\mathbf{c}^n))} \le 2 \|\nabla \mathbf{u}\|_{L^\infty} \|\mathbf{c}^{n+1}\|_{L^1}, $$

so, thanks to Lemma B.1,

$$ \frac{1}{\delta t_0} (b - \operatorname{tr}(\mathbf{c}^n)) + \frac{b(b-d)}{\lambda(b - \operatorname{tr}(\mathbf{c}^n))} \le 2d \|\nabla \mathbf{u}\|_{L^\infty} \operatorname{tr}(\mathbf{c}^{n+1}) = 2db \|\nabla \mathbf{u}\|_{L^\infty}. \quad (B.11) $$

We observe that the two terms at the left-hand side are positive. This inequality is thus impossible, and we have $\operatorname{tr}(\mathbf{c}^{n+1}) < b$ for any value of the time step, if

$$ \frac{b-d}{2d\lambda(b - \operatorname{tr}(\mathbf{c}^n))} > \|\nabla \mathbf{u}\|_{L^\infty}. \quad (B.12) $$
Let us suppose that this condition is not satisfied. Inequality (B.11) also yields:

$$ \frac{1}{\delta t_0} (b - \text{tr}(c^n)) \le 2 \, db \, ||\nabla u||_{L^\infty}, $$

which, by contradiction, shows that $\text{tr}(c^{n+1}) < b$ provided that

$$ \delta t < \frac{b - \text{tr}(c^n)}{2 \, db \, ||\nabla u||_{L^\infty}}. \qquad (B.13) $$

Finally, keeping all the terms of the left-hand side of (B.11), the same argument yields that $\text{tr}(c^{n+1}) < b$ if

$$ \delta t < \frac{b - \text{tr}(c^n)}{2 \, db \, ||\nabla u||_{L^\infty} - \frac{b(b-d)}{\lambda(b - \text{tr}(c^n))}}. $$

The right-hand side may be seen as a function of the quantity $b - \text{tr}(c^n)$ defined on the part of $\mathbb{R}$ where the denominator is positive, and which takes a minimum value. This latter thus gives an upper bound for $\delta t$ which is independent of $c^n$ and reads:

$$ \delta t < \frac{b-d}{d^2 b \lambda ||\nabla u||_{L^\infty}^2}. \qquad (B.14) $$

Summing up the developments of this section, we get the following result.

**Lemma B.2.** Let $\delta t$ satisfy (B.10), and $c^n$ be a symmetric definite positive tensor such that $\text{tr}(c^n) < b$. Then System (B.9) has a unique solution $c^{n+1}$, which is a symmetric definite positive tensor. If, in addition, one of the inequalities (B.12), (B.13) or (B.14) is satisfied, then $\text{tr}(c^{n+1}) < b$.

## APPENDIX C. A STABILITY PROPERTY OF THE UPWIND CONVECTION OF TENSORS

Let $c \in \mathbb{R}^{d \times d}$, $u$ be a divergence-free regular velocity field and let $\mathcal{F}: \mathbb{R}^{d \times d} \to \mathbb{R}$ be a regular function. Then, we have

$$ (\mathcal{F}'(c), \partial_t c + \text{div}(u c)) = \partial_t \mathcal{F}(c) + \text{div}(u \mathcal{F}(c)). $$

The aim of this section is to state and prove a discrete analogue of this relation with an implicit upwind discretization of the convection operator. This result turns out to be an inequality and holds only when the function $\mathcal{F}$ is convex, as stated in Lemma C.1 below. A similar result for a scalar unknown is an easy consequence of Lemma A2 from [30]. The specialization of Lemma C.1 to the function of interest in the stability analysis of the scheme, namely $c \mapsto -\text{tr}(\ln c)$, is used and proved in [14].

**Lemma C.1.** Let $\mathcal{C}$ be an open convex subset of $\mathbb{R}^{d \times d}$, and let $\mathcal{F}$ a convex differentiable function from $\mathcal{C}$ to $\mathbb{R}$. Let $c$ and $c^*$ be two elements of $\mathcal{C}$, $\alpha$ a positive real number and $(u_\sigma)_{\sigma \in \mathcal{E}}$ a family of real numbers ($\mathcal{E}$ is a finite set of indexes) satisfying

$$ \sum_{\sigma \in \mathcal{E}} u_{\sigma} = 0. \qquad (C.15) $$

For $\sigma \in \mathcal{E}$, let $c_\sigma \in \mathcal{C}$ and let us suppose that $c_\sigma = c$ if $u_\sigma \ge 0$. Then

$$ \left( \mathcal{F}'(c), \alpha(c - c^*) + \sum_{\sigma \in \mathcal{E}} u_{\sigma} c_{\sigma} \right) \ge \alpha(\mathcal{F}(c) - \mathcal{F}(c^*)) + \sum_{\sigma \in \mathcal{E}} u_{\sigma} \mathcal{F}(c_{\sigma}). $$
Let us suppose that this condition is not satisfied. Inequality (B.11) also yields:

$$ \frac{1}{\delta t_0} (b - \text{tr}(c^n)) \le 2 \, db \, ||\nabla u||_{L^\infty}, $$

which, by contradiction, shows that $\text{tr}(c^{n+1}) < b$ provided that

$$ \delta t < \frac{b - \text{tr}(c^n)}{2 \, db \, ||\nabla u||_{L^\infty}}. \qquad (B.13) $$

Finally, keeping all the terms of the left-hand side of (B.11), the same argument yields that $\text{tr}(c^{n+1}) < b$ if

$$ \delta t < \frac{b - \text{tr}(c^n)}{2 \, db \, ||\nabla u||_{L^\infty} - \frac{b(b-d)}{\lambda(b - \text{tr}(c^n))}}. $$

The right-hand side may be seen as a function of the quantity $b - \text{tr}(c^n)$ defined on the part of $\mathbb{R}$ where the denominator is positive, and which takes a minimum value. This latter thus gives an upper bound for $\delta t$ which is independent of $c^n$ and reads:

$$ \delta t < \frac{b-d}{d^2 b \lambda ||\nabla u||_{L^\infty}^2}. \qquad (B.14) $$

Summing up the developments of this section, we get the following result.

**Lemma B.2.** Let $\delta t$ satisfy (B.10), and $c^n$ be a symmetric definite positive tensor such that $\text{tr}(c^n) < b$. Then System (B.9) has a unique solution $c^{n+1}$, which is a symmetric definite positive tensor. If, in addition, one of the inequalities (B.12), (B.13) or (B.14) is satisfied, then $\text{tr}(c^{n+1}) < b$.

## APPENDIX C. A STABILITY PROPERTY OF THE UPWIND CONVECTION OF TENSORS

Let $c \in \mathbb{R}^{d \times d}$, $u$ be a divergence-free regular velocity field and let $\mathcal{F}: \mathbb{R}^{d \times d} \to \mathbb{R}$ be a regular function. Then, we have

$$ (\mathcal{F}'(c), \partial_t c + \text{div}(u c)) = \partial_t \mathcal{F}(c) + \text{div}(u \mathcal{F}(c)). $$

The aim of this section is to state and prove a discrete analogue of this relation with an implicit upwind discretization of the convection operator. This result turns out to be an inequality and holds only when the function $\mathcal{F}$ is convex, as stated in Lemma C.1 below. A similar result for a scalar unknown is an easy consequence of Lemma A2 from [30]. The specialization of Lemma C.1 to the function of interest in the stability analysis of the scheme, namely $c \mapsto -\text{tr}(\ln c)$, is used and proved in [14].

**Lemma C.1.** Let $\mathcal{C}$ be an open convex subset of $\mathbb{R}^{d \times d}$, and let $\mathcal{F}$ a convex differentiable function from $\mathcal{C}$ to $\mathbb{R}$. Let $c$ and $c^*$ be two elements of $\mathcal{C}$, $\alpha$ a positive real number and $(u_\sigma)_{\sigma \in \mathcal{E}}$ a family of real numbers ($\mathcal{E}$ is a finite set of indexes) satisfying

$$ \sum_{\sigma \in \mathcal{E}} u_{\sigma} = 0. \qquad (C.15) $$

For $\sigma \in \mathcal{E}$, let $c_\sigma \in \mathcal{C}$ and let us suppose that $c_\sigma = c$ if $u_\sigma \ge 0$. Then

$$ \left( \mathcal{F}'(c), \alpha(c - c^*) + \sum_{\sigma \in \mathcal{E}} u_{\sigma} c_{\sigma} \right) \ge \alpha(\mathcal{F}(c) - \mathcal{F}(c^*)) + \sum_{\sigma \in \mathcal{E}} u_{\sigma} \mathcal{F}(c_{\sigma}). $$
Proof. Let us write $(\mathcal{F}'(c), \alpha(c - c^*) + \sum_{\sigma \in \mathcal{E}} u_{\sigma}c_{\sigma}) = T_{\partial/\partial_t} + T_{\text{div}}$ with

$$T_{\partial/\partial_t} = (\mathcal{F}'(c), \alpha(c - c^*)) \quad \text{and} \quad T_{\text{div}} = \left( \mathcal{F}'(c), \sum_{\sigma \in \mathcal{E}} u_{\sigma}c_{\sigma} \right).$$

By convexity of the function $\mathcal{F}$, we have $T_{\partial/\partial_t} \ge \alpha(\mathcal{F}(c) - \mathcal{F}(c^*))$. In addition, thanks to the equality (C.15), we may recast $T_{\text{div}}$ as

$$T_{\text{div}} = \sum_{\sigma \in \mathcal{E}} u_{\sigma} \mathcal{F}(c_{\sigma}) + \sum_{\sigma \in \mathcal{E}} u_{\sigma} R_{\sigma} \quad \text{with} \quad R_{\sigma} = (\mathcal{F}'(c), c_{\sigma} - c) + \mathcal{F}(c) - \mathcal{F}(c_{\sigma}).$$

By convexity of $\mathcal{F}$, $R_{\sigma} \le 0$. In addition, when $u_{\sigma} \ge 0$, we have $c_{\sigma} = c$ (by assumption) and $R_{\sigma}$ vanishes. The product $u_{\sigma} R_{\sigma}$ is thus non-negative, which concludes the proof. $\square$

Let us now apply this lemma to provide the argument requested in Section 3.2 for the stability of the scheme, i.e. a discrete analogue of relation (19). First, we recall the following classical results for functions of tensors [15].

**Lemma C.2.** Let $\mathcal{F}$ be the function $S_{+}^*(\mathbb{R}^{d \times d}) \to \mathbb{R}$ defined by $c \mapsto -\ln \det(c)$. Then $\mathcal{F}$ is convex, regular and its first derivative is

$$(\mathcal{F}'(c), d) = -\operatorname{tr}(c^{-1} d), \quad \forall d \in S_{+}^*(\mathbb{R}^{d \times d}).$$

In addition, for any $c \in S_{+}^*(\mathbb{R}^{d \times d})$, $\operatorname{tr}(\ln c) = \ln \det(c)$.

Returning to the scheme for the transport of $c$, the application of Lemma C.1 to the function $\mathcal{F}(c) = -\operatorname{tr}(\ln c)$ readily yields

$$\forall K \in \mathcal{M}, \quad \operatorname{tr} \left( (c_K^{n+1})^{-1} \left( \frac{|K|}{\delta t} (c_K^{n+1} - c_K^n) + \operatorname{div}_K (u^n c^{n+1}) \right) \right) \\ \le \frac{|K|}{\delta t} \left( \operatorname{tr}(\ln c_K^{n+1}) - \operatorname{tr}(\ln c_K^n) \right) + \operatorname{div}_K (u^n \operatorname{tr}(\ln c^{n+1})).$$

This is the identity needed to prove the stability of the scheme.

## APPENDIX D. A SCALAR MODEL PROBLEM

We consider the following model problem on $I \times (0, T)$, with $I = (-1, 1)$,

$$\partial_t c + \partial_x c = a(x) c \text{ on } I \times (0, T), \quad c(x, 0) = 1 \text{ for } x \in I, \quad c(-1, t) = c_L \text{ for } t \in (0, T). \qquad (\text{D.16})$$

Since the inlet condition $c_L$ is constant, the solution of this problem is stationary for $t \ge 2$. The form of the source term $a(x)c$ is chosen to mimick the production term $(\nabla u)c + c(\nabla u)^t$ involved in the transport equation of the conformation tensor. In addition, we are interested here in the situation where $a(x)$ takes locally large values, so as to simulate zones of high elongation in a viscoelastic flows. Our understanding of viscoelastic flows is that the stiffness of this right-hand side is actually the origin of the so-called High Weissenberg Problem, i.e. the difficulty for numerical schemes to accurately capture (or, worse, not to blow up) for moderate-to-high Weissenberg numbers. Finally, even if it is probably not essential, we chose $a(x)$ in such a way that $c(1, t) = c_L$ which, to some extent, may be considered to represent what happens along a closed streamline in a steady viscoelastic flow. The chosen source term $a(x)$ and the corresponding stationary solution $c_s(x)$ are

$$a(x) = \begin{cases} 0 & -1 \le x < -\delta, \\ -\beta \sin(\frac{\pi}{\delta} x) & -\delta \le x \le \delta, \\ 0 & \delta < x \le 1. \end{cases} \quad c_s(x) = \begin{cases} c_L & -1 \le x < -\delta, \\ c_L \exp\left(\beta \frac{\delta}{\pi} \left(1 + \cos(\frac{\pi}{\delta} x)\right)\right) & -\delta \le x \le \delta, \\ c_L & \delta < x \le 1. \end{cases}$$
Proof. Let us write $(\mathcal{F}'(c), \alpha(c - c^*) + \sum_{\sigma \in \mathcal{E}} u_{\sigma} c_{\sigma}) = T_{\partial/\partial_t} + T_{\text{div}}$ with

$$T_{\partial/\partial_t} = (\mathcal{F}'(c), \alpha(c - c^*)) \quad \text{and} \quad T_{\text{div}} = \left( \mathcal{F}'(c), \sum_{\sigma \in \mathcal{E}} u_{\sigma} c_{\sigma} \right).$$

By convexity of the function $\mathcal{F}$, we have $T_{\partial/\partial_t} \ge \alpha(\mathcal{F}(c) - \mathcal{F}(c^*))$. In addition, thanks to the equality (C.15), we may recast $T_{\text{div}}$ as

$$T_{\text{div}} = \sum_{\sigma \in \mathcal{E}} u_{\sigma} \mathcal{F}(c_{\sigma}) + \sum_{\sigma \in \mathcal{E}} u_{\sigma} R_{\sigma} \quad \text{with} \quad R_{\sigma} = (\mathcal{F}'(c), c_{\sigma} - c) + \mathcal{F}(c) - \mathcal{F}(c_{\sigma}).$$

By convexity of $\mathcal{F}$, $R_{\sigma} \le 0$. In addition, when $u_{\sigma} \ge 0$, we have $c_{\sigma} = c$ (by assumption) and $R_{\sigma}$ vanishes. The product $u_{\sigma} R_{\sigma}$ is thus non-negative, which concludes the proof. $\square$

Let us now apply this lemma to provide the argument requested in Section 3.2 for the stability of the scheme, i.e. a discrete analogue of relation (19). First, we recall the following classical results for functions of tensors [15].

**Lemma C.2.** Let $\mathcal{F}$ be the function $\mathbb{S}_+^*(\mathbb{R}^{d \times d}) \to \mathbb{R}$ defined by $c \mapsto -\ln \det(c)$. Then $\mathcal{F}$ is convex, regular and its first derivative is

$$(\mathcal{F}'(c), d) = -\operatorname{tr}(c^{-1} d), \quad \forall d \in \mathbb{S}_+^*(\mathbb{R}^{d \times d}).$$

In addition, for any $c \in \mathbb{S}_+^*(\mathbb{R}^{d \times d})$, $\operatorname{tr}(\ln c) = \ln \det(c)$.

Returning to the scheme for the transport of $c$, the application of Lemma C.1 to the function $\mathcal{F}(c) = -\operatorname{tr}(\ln c)$ readily yields

$$\forall K \in \mathcal{M}, \quad \operatorname{tr} \left( (c_K^{n+1})^{-1} \left( \frac{|K|}{\delta t} (c_K^{n+1} - c_K^n) + \operatorname{div}_K (u^n c^{n+1}) \right) \right) \\ \le \frac{|K|}{\delta t} \left( \operatorname{tr}(\ln c_K^{n+1}) - \operatorname{tr}(\ln c_K^n) \right) + \operatorname{div}_K (u^n \operatorname{tr}(\ln c^{n+1})).$$

This is the identity needed to prove the stability of the scheme.

## APPENDIX D. A SCALAR MODEL PROBLEM

We consider the following model problem on $I \times (0, T)$, with $I = (-1, 1)$,

$$\partial_t c + \partial_x c = a(x) c \text{ on } I \times (0, T), \quad c(x, 0) = 1 \text{ for } x \in I, \quad c(-1, t) = c_L \text{ for } t \in (0, T). \qquad (\text{D.16})$$

Since the inlet condition $c_L$ is constant, the solution of this problem is stationary for $t \ge 2$. The form of the source term $a(x)c$ is chosen to mimick the production term $(\nabla u)c + c(\nabla u)^t$ involved in the transport equation of the conformation tensor. In addition, we are interested here in the situation where $a(x)$ takes locally large values, so as to simulate zones of high elongation in a viscoelastic flows. Our understanding of viscoelastic flows is that the stiffness of this right-hand side is actually the origin of the so-called High Weissenberg Problem, i.e. the difficulty for numerical schemes to accurately capture (or, worse, not to blow up) for moderate-to-high Weissenberg numbers. Finally, even if it is probably not essential, we chose $a(x)$ in such a way that $c(1, t) = c_L$ which, to some extent, may be considered to represent what happens along a closed streamline in a steady viscoelastic flow. The chosen source term $a(x)$ and the corresponding stationary solution $c_s(x)$ are

$$a(x) = \begin{cases} 0 & -1 \le x < -\delta, \\ -\beta \sin(\frac{\pi}{\delta} x) & -\delta \le x \le \delta, \\ 0 & \delta < x \le 1. \end{cases} \qquad c_s(x) = \begin{cases} c_L & -1 \le x < -\delta, \\ c_L \exp\left(\beta \frac{\delta}{\pi} \left(1 + \cos(\frac{\pi}{\delta} x)\right)\right) & -\delta \le x \le \delta, \\ c_L & \delta < x \le 1. \end{cases}$$
where $\beta$ and $\delta < 1$ are given parameters. The problem may be made stiffer by decreasing $\delta$ and increasing $\beta$.

Let us consider a uniform partition $0 = t_0 < t_1 < \dots < t_N = T$ of the time interval $(0, T)$, and let $\delta t = t_{n+1} - t_n$ for $0 \le n \le N - 1$ be the constant time step. Let us also consider the uniform partition $-1 = x_0 < x_1 < \dots < x_K = 1$ of the space interval $I$ with $h = x_{k+1} - x_k$ for $0 \le k \le K - 1$ the constant space step. The solution of (D.16) is approximated by piecewise constant unknowns $c_{k+\frac{1}{2}}^n$, for $0 \le k \le K - 1$ and $0 \le n \le N$, where $c_{k+\frac{1}{2}}^n$ is an approximation of the mean value of $c(x, t_n)$ over $(x_k, x_{k+1})$. The scheme is initialized by

$$c_{k+\frac{1}{2}}^0 = 1, \quad \text{for } 0 \le k \le K-1.$$

Then the unknowns are updated in time by an operator-splitting technique, which reads, in a semi-discrete setting

**Advection I – Solve for $c^{n+\frac{1}{3}}$:**

$$\partial_t c + \partial_x c = 0 \text{ on } I \times \left(t_n, t_n + \frac{\delta t}{2}\right), \quad c(x, t_n) = c^n \text{ for } x \in I. \qquad (D.17a)$$

**ODE – Set $c_{k+\frac{1}{2}}(t_n) = c_{k+\frac{1}{2}}^{n+\frac{1}{3}}$ and solve for $c_{k+\frac{1}{2}}^{n+\frac{2}{3}} = c_{k+\frac{1}{2}}(t_n + \delta t)$:**

$$\partial_t c_{k+\frac{1}{2}}(t) = a_{k+\frac{1}{2}} c_{k+\frac{1}{2}}(t), \quad \text{for } 0 \le k \le K-1. \qquad (D.17b)$$

**Advection II – Solve for $c^{n+1}$:**

$$\partial_t c + \partial_x c = 0 \text{ on } I \times \left(t_n + \frac{\delta t}{2}, t_{n+1}\right), \quad c\left(x, t_n + \frac{\delta t}{2}\right) = c^{n+\frac{2}{3}} \text{ for } x \in I. \qquad (D.17c)$$

In the ODE step, the real number $a_{k+\frac{1}{2}}$ is defined as the mean value of $a(x)$ over $(x_k, x_{k+1})$. This step is solved quasi-exactly by

$$c_{k+\frac{1}{2}}^{n+\frac{2}{3}} = c_{k+\frac{1}{2}}^{n+\frac{1}{3}} \exp(a_{k+\frac{1}{2}} \delta t).$$

We test four different ways of solving the advection steps:

(i) The usual implicit first order upwind scheme

$$c_{k+\frac{1}{2}}^{n+1} = c_{k+\frac{1}{2}}^{n} - \frac{\delta t}{h} \left( c_{k+\frac{1}{2}}^{n+1} - c_{k}^{n+1} \right), \quad \text{for } 0 \le k \le K-1,$$

where $c_k^{n+1}$ stands for $c_{k-\frac{1}{2}}^{n+1}$ for $k \ge 1$ and $c_0^{n+1} = c_L$.

(ii) The same upwind scheme applied to the following logarithmic formulation of the transport equation (easily obtained by multiplication of equation (D.16) by $1/c$)

$$\log(c_{k+\frac{1}{2}}^{n+1}) = \log(c_{k+\frac{1}{2}}^{n}) - \frac{\delta t}{h} \left( \log(c_{k+\frac{1}{2}}^{n+1}) - \log(c_{k}^{n+1}) \right), \quad \text{for } 0 \le k \le K-1.$$

(iii) An explicit MUSCL-type scheme which reads

$$c_{k+\frac{1}{2}}^{n+1} = c_{k+\frac{1}{2}}^{n} - \frac{\delta t}{h} (c_{k+1}^{n} - c_{k}^{n}), \quad \text{for } 0 \le k \le K-1,$$

where $c_{k+1}^n$, the approximation of the unknown at $x_{k+1}^n$, is defined, for $1 \le k \le K-2$ by

$$c_{k+1}^{n} - c_{k+\frac{1}{2}}^{n} = \min\min\left(\left(c_{k+\frac{1}{2}}^{n} - c_{k+\frac{1}{2}}^{n}\right)/2, c_{k+\frac{1}{2}}^{n} - c_{k-\frac{1}{2}}^{n}\right),$$

with the usual definition of the minmod function, namely, for any pair of real numbers $(a, b)$, $\min\min(a, b) = \operatorname{sgn}(a, b) \min(|a|, |b|)$, $\operatorname{sgn}(a, b)$ being equal to the sign of $a$ and $b$ if $a$ and $b$ have the same sign and 0
where $\beta$ and $\delta < 1$ are given parameters. The problem may be made stiffer by decreasing $\delta$ and increasing $\beta$.

Let us consider a uniform partition $0 = t_0 < t_1 < \cdots < t_N = T$ of the time interval $(0, T)$, and let $\delta t = t_{n+1} - t_n$ for $0 \le n \le N - 1$ be the constant time step. Let us also consider the uniform partition $-1 = x_0 < x_1 < \cdots < x_K = 1$ of the space interval $I$ with $h = x_{k+1} - x_k$ for $0 \le k \le K - 1$ the constant space step. The solution of (D.16) is approximated by piecewise constant unknowns $c_{k+\frac{1}{2}}^n$, for $0 \le k \le K - 1$ and $0 \le n \le N$, where $c_{k+\frac{1}{2}}^n$ is an approximation of the mean value of $c(x, t_n)$ over $(x_k, x_{k+1})$. The scheme is initialized by

$$c_{k+\frac{1}{2}}^0 = 1, \quad \text{for } 0 \le k \le K-1.$$

Then the unknowns are updated in time by an operator-splitting technique, which reads, in a semi-discrete setting

Advection I – Solve for $c^{n+\frac{1}{3}}$:

$$\partial_t c + \partial_x c = 0 \text{ on } I \times \left(t_n, t_n + \frac{\delta t}{2}\right), \quad c(x, t_n) = c^n \text{ for } x \in I. \qquad (D.17a)$$

ODE – Set $c_{k+\frac{1}{2}}(t_n) = c_{k+\frac{1}{2}}^{n+\frac{1}{3}}$ and solve for $c_{k+\frac{1}{2}}^{n+\frac{2}{3}} = c_{k+\frac{1}{2}}(t_n + \delta t)$:

$$\partial_t c_{k+\frac{1}{2}}(t) = a_{k+\frac{1}{2}} c_{k+\frac{1}{2}}(t), \quad \text{for } 0 \le k \le K-1. \qquad (D.17b)$$

Advection II – Solve for $c^{n+1}$:

$$\partial_t c + \partial_x c = 0 \text{ on } I \times \left(t_n + \frac{\delta t}{2}, t_{n+1}\right), \quad c\left(x, t_n + \frac{\delta t}{2}\right) = c^{n+\frac{2}{3}} \text{ for } x \in I. \qquad (D.17c)$$

In the ODE step, the real number $a_{k+\frac{1}{2}}$ is defined as the mean value of $a(x)$ over $(x_k, x_{k+1})$. This step is solved quasi-exactly by

$$c_{k+\frac{1}{2}}^{n+\frac{2}{3}} = c_{k+\frac{1}{2}}^{n+\frac{1}{3}} \exp(a_{k+\frac{1}{2}} \delta t).$$

We test four different ways of solving the advection steps:

(i) The usual implicit first order upwind scheme

$$c_{k+\frac{1}{2}}^{n+1} = c_{k+\frac{1}{2}}^{n} - \frac{\delta t}{h} \left( c_{k+\frac{1}{2}}^{n+1} - c_{k}^{n+1} \right), \quad \text{for } 0 \le k \le K-1,$$

where $c_k^{n+1}$ stands for $c_{k-\frac{1}{2}}^{n+1}$ for $k \ge 1$ and $c_0^{n+1} = c_L$.

(ii) The same upwind scheme applied to the following logarithmic formulation of the transport equation (easily obtained by multiplication of equation (D.16) by $1/c$)

$$\log(c_{k+\frac{1}{2}}^{n+1}) = \log(c_{k+\frac{1}{2}}^{n}) - \frac{\delta t}{h} \left( \log(c_{k+\frac{1}{2}}^{n+1}) - \log(c_{k}^{n+1}) \right), \quad \text{for } 0 \le k \le K-1.$$

(iii) An explicit MUSCL-type scheme which reads

$$c_{k+\frac{1}{2}}^{n+1} = c_{k+\frac{1}{2}}^{n} - \frac{\delta t}{h} (c_{k+1}^{n} - c_{k}^{n}), \quad \text{for } 0 \le k \le K-1,$$

where $c_{k+1}^n$, the approximation of the unknown at $x_{k+1}^n$, is defined, for $1 \le k \le K-2$ by

$$c_{k+1}^{n} - c_{k+\frac{1}{2}}^{n} = \minmod\left(\left(c_{k+\frac{3}{2}}^{n} - c_{k+\frac{1}{2}}^{n}\right)/2, c_{k+\frac{1}{2}}^{n} - c_{k-\frac{1}{2}}^{n}\right),$$

with the usual definition of the minmod function, namely, for any pair of real numbers $(a, b)$, $\minmod(a, b) = \operatorname{sgn}(a, b) \min(|a|, |b|)$, $\operatorname{sgn}(a, b)$ being equal to the sign of $a$ and $b$ if $a$ and $b$ have the same sign and 0
<figure><img src="image_31.png" /><figcaption>FIGURE D.1. Comparison of the different methods of discretization of the advection step with $\beta = 50$ and $\delta = 0.1$. Left: stationary states obtained for $K = 200$. Right: relative errors (in discrete $L^2$-norm) with respect to the analytical solution as a function of the number of cells.</figcaption></figure>

otherwise. The definition of the face approximations for $k = 0$, $k = 1$ and $k = K$ is the same as for the upwind scheme (which, in fact, does not matter in the application case below, since the solution is almost constant near the beginning and end of the interval $I$).

(iv) The same MUSCL-type scheme applied to the logarithmic formulation of choice (ii).

We compare the resolution of equation (D.16) with these different numerical schemes. The stiffness parameters are fixed to $\beta = 50$ and $\delta = 0.1$, which indeed yields a rather stiff problem. The time step is fixed to $\delta t = 1/(2K)$, so as to satisfy the stability constraints of the (explicit-in-time) MUSCL schemes. Results for $K = 200$ are plotted in Figure D.1, together with a convergence study with respect to the number of cells. In both cases, we consider the results at $t = 10$, when the solution is stationary both in the continuous and discrete cases. The first conclusion is that, even with the accurate solution of the ODE step implemented here, large errors may be obtained. Second, we observe that switching to the logarithmic formulation in the transport steps dramatically improves the scheme accuracy. On the opposite, the MUSCL discretization does not bring a decisive gain in accuracy. These conclusions motivate our choice for the scheme developed in this paper for the transport of the conformation tensor: we use the logarithmic formulation for the transport steps and keep an implicit linear and monotone (so first order upwind) scheme for their solution.

Note that, for one-dimensional problems as considered in this section, more sophisticated methods to cope with stiff source terms are proposed in the literature such as the so-called “well-balanced” schemes, designed to preserve at the discrete level specific “equilibrium” (i.e. steady) states [11, 26, 33]. Such schemes are not tested here because their extension to multi-dimensional flow problems is not straightforward.

*Acknowledgements.* The authors would like to thank R. de Loubens, from TotalEnergies, for valuable discussions on this topic, F. Babik, from the CALIF³S development team at IRSN, for his help in implementing the scheme and TotalEnergies for funding this work and for providing access to the Pangea II supercomputer.

## REFERENCES

[1] A.M. Afonso, F.T. Pinho and M.A. Alves, The kernel-conformation constitutive laws. *J. Non-Newtonian Fluid Mech.* **167** (2012) 30–37.

[2] M.A. Alves, P.J. Oliveira and F.T. Pinho, Numerical methods for viscoelastic fluid flows. *Ann. Rev. Fluid Mech.* **53** (2021) 509–541.

[3] E. Audusse, F. Bouchut, M.-O. Bristeau, R. Klein and B. Perthame, A fast and stable well-balanced scheme with hydrostatic reconstruction for Shallow Water flows. *SIAM J. Sci. Comput.* **25** (2004) 2050–2065.

[4] N. Balci, B. Thomases, M. Renardy and C.R. Doering, Symmetric factorization of the conformation tensor in viscoelastic fluid models. *J. Non-Newtonian Fluid Mech.* **166** (2011) 546–553.
<figure><img src="image_31.png" /><figcaption>FIGURE D.1. Comparison of the different methods of discretization of the advection step with $\beta = 50$ and $\delta = 0.1$. Left: stationary states obtained for $K = 200$. Right: relative errors (in discrete $L^2$-norm) with respect to the analytical solution as a function of the number of cells.</figcaption></figure>

otherwise. The definition of the face approximations for $k = 0$, $k = 1$ and $k = K$ is the same as for the upwind scheme (which, in fact, does not matter in the application case below, since the solution is almost constant near the beginning and end of the interval $I$).

(iv) The same MUSCL-type scheme applied to the logarithmic formulation of choice (ii).

We compare the resolution of equation (D.16) with these different numerical schemes. The stiffness parameters are fixed to $\beta = 50$ and $\delta = 0.1$, which indeed yields a rather stiff problem. The time step is fixed to $\delta t = 1/(2K)$, so as to satisfy the stability constraints of the (explicit-in-time) MUSCL schemes. Results for $K = 200$ are plotted in Figure D.1, together with a convergence study with respect to the number of cells. In both cases, we consider the results at $t = 10$, when the solution is stationary both in the continuous and discrete cases. The first conclusion is that, even with the accurate solution of the ODE step implemented here, large errors may be obtained. Second, we observe that switching to the logarithmic formulation in the transport steps dramatically improves the scheme accuracy. On the opposite, the MUSCL discretization does not bring a decisive gain in accuracy. These conclusions motivate our choice for the scheme developed in this paper for the transport of the conformation tensor: we use the logarithmic formulation for the transport steps and keep an implicit linear and monotone (so first order upwind) scheme for their solution.

Note that, for one-dimensional problems as considered in this section, more sophisticated methods to cope with stiff source terms are proposed in the literature such as the so-called “well-balanced” schemes, designed to preserve at the discrete level specific “equilibrium” (i.e. steady) states [11, 26, 33]. Such schemes are not tested here because their extension to multi-dimensional flow problems is not straightforward.

*Acknowledgements.* The authors would like to thank R. de Loubens, from TotalEnergies, for valuable discussions on this topic, F. Babik, from the CALIF³S development team at IRSN, for his help in implementing the scheme and TotalEnergies for funding this work and for providing access to the Pangea II supercomputer.

## REFERENCES

[1] A.M. Afonso, F.T. Pinho and M.A. Alves, The kernel-conformation constitutive laws. *J. Non-Newtonian Fluid Mech.* **167** (2012) 30–37.

[2] M.A. Alves, P.J. Oliveira and F.T. Pinho, Numerical methods for viscoelastic fluid flows. *Ann. Rev. Fluid Mech.* **53** (2021) 509–541.

[3] E. Audusse, F. Bouchut, M.-O. Bristeau, R. Klein and B. Perthame, A fast and stable well-balanced scheme with hydrostatic reconstruction for Shallow Water flows. *SIAM J. Sci. Comput.* **25** (2004) 2050–2065.

[4] N. Balci, B. Thomases, M. Renardy and C.R. Doering, Symmetric factorization of the conformation tensor in viscoelastic fluid models. *J. Non-Newtonian Fluid Mech.* **166** (2011) 546–553.
[5] J.W. Barrett and S. Boyaval, Existence and approximation of a (regularized) Oldroyd-B model. *Math. Models Methods Appl. Sci.* **21** (2011) 1783–1837.

[6] J.W. Barrett and S. Boyaval, Finite element approximation of the FENE-P model. *IMA J. Numer. Anal.* **38** (2018) 1599–1660.

[7] J.W. Barrett and E. Süli, Existence of global weak solutions to some regularized kinetic models for dilute polymers. *Multiscale Model. Simul.* **6** (2007) 506–546.

[8] J.W. Barrett, C. Schwab and E. Süli, Existence of global weak solutions for some polymeric flow models. *Math. Models Methods Appl. Sci.* **15** (2005) 939–983.

[9] A. Bermudez and M.E. Vazquez, Upwind methods for hyperbolic conservation laws with source terms. *Comput. Fluids* **23** (1994) 1049–1071.

[10] R.B. Bird and J.M. Wiest, Constitutive equations for polymeric liquids. *Ann. Rev. Fluid Mech.* **27** (1995) 169–193.

[11] R. Botchorishvili, B. Perthame and A. Vasseur, Equilibrium schemes for scalar conservation laws with stiff sources. *Math. Comput.* **72** (2003) 131–157.

[12] F. Bouchut, H. Ounaissa and B. Perthame, Upwinding of the source term at interfaces for Euler equations with high friction. *Comput. Math. App.* **53** (2007) 361–375.

[13] S. Boyaval, Lid-driven-cavity simulations of Oldroyd-B models using free-energy-dissipative schemes, in Numerical Mathematics and Advanced Applications 2009. Springer (2010) 191–198.

[14] S. Boyaval, T. Lelièvre and C. Mangoubi, Free-energy-dissipative schemes for the Oldroyd-B model. *Math. Model. Numer. Anal.* **43** (2009) 523–561.

[15] S. Boyd and L. Vandenberghe, Convex Optimization. Cambridge University Press, Cambridge (2004).

[16] CALIF³S, A software components library for the computation of fluid flows. https://gforge.irsn.fr/gf/project/califs.

[17] G.-Q. Chen, C. D. Levermore and T.-P. Liu, Hyperbolic conservation laws with stiff relaxation terms and entropy. *Commun. Pure Appl. Math.* **47** (1994) 787–830.

[18] R. Comminal, J. Spangenberg and J.H. Hattel, Robust simulations of viscoelastic flows at high Weissenberg numbers with the streamfunction/log-conformation formulation. *J. Non-Newtonian Fluid Mech.* **223** (2015) 37–61.

[19] M. Crouzeix and P.A. Raviart, Conforming and nonconforming finite element methods for solving the stationary Stokes equations. *RAIRO Série Rouge* **7** (1973) 33–75.

[20] H. Damanik, J. Hron, A. Ouazzi and S. Turek, A monolithic FEM approach for the log-conformation reformulation (LCR) of viscoelastic flow problems. *J. Non-Newtonian Fluid Mech.* **165** (2010) 1105–1113.

[21] A. Duran, F. Marche, R. Turpault and C. Berthon, Asymptotic preserving scheme for the shallow water equations with source terms on unstructured meshes. *J. Comput. Phys.* **287** (2015) 184–206.

[22] R. Fattal and R. Kupferman, Constitutive laws for the matrix-logarithm of the conformation tensor. *J. Non-Newtonian Fluid Mech.* **123** (2004) 281–285.

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

[24] T. Gallouët, R. Herbin, J.-C. Latché and K. Mallem, Convergence of the Marker-And-Cell scheme for the incompressible Navier-Stokes equations on non-uniform grids. *Found. Comput. Math.* **18** (2018) 249–289.

[25] D. Grapsas, R. Herbin, W. Kheriji and J.-C. Latché, An unconditionally stable staggered pressure correction scheme for the compressible Navier-Stokes equations. *SMAI J. Comput. Math.* **2** (2016) 51–97.

[26] J.M. Greenberg and A.Y. Leroux, A well-balanced scheme for the numerical processing of source terms in hyperbolic equations. *SIAM J. Numer. Anal.* **33** (1996) 1–16.

[27] A.M. Grillet, B. Yang, B. Khomami and E.S. Shaqfeh, Modeling of viscoelastic lid driven cavity flow using finite element simulations. *J. Non-Newtonian Fluid Mech.* **88** (1999) 99–131.

[28] F.H. Harlow and A.A. Amsden, A numerical fluid dynamics calculation method for all flow speeds. *J. Comput. Phys.* **8** (1971) 197–213.

[29] F.H. Harlow and J.E. Welsh, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface. *Phys. Fluids* **8** (1965) 2182–2189.

[30] R. Herbin, W. Kheriji and J.-C. Latché, On some implicit and semi-implicit staggered schemes for the shallow water and Euler equations. *Math. Model. Numer. Anal.* **48** (2014) 1807–1857.

[31] D. Hu and T. Lelievre, New entropy estimates for the Oldroyd-B model and related models. *Commun. Math. Sci.* **5** (2007) 909–916.

[32] M.A. Hulsen, R. Fattal and R. Kupferman, Flow of viscoelastic fluids past a cylinder at high Weissenberg number: stabilized simulations using matrix logarithms. *J. Non-Newtonian Fluid Mech.* **127** (2005) 27–39.

[33] S. Jin, Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations. *SIAM J. Sci. Comput.* **21** (1999) 441–454.

[34] T. Katsaounis, B. Perthame and C. Simeoni, Upwinding sources at interfaces in conservation laws. *Appl. Math. Lett.* **17** (2004) 309–316.

[35] J.-C. Latché and K. Saleh, A convergent staggered scheme for the variable density incompressible Navier-Stokes equations. *Math. Comput.* **87** (2018) 581–632.

[36] J.-C. Latché, B. Piar and K. Saleh, A discrete kinetic energy preserving convection operator for variable density flows on locally refined staggered meshes. https://hal.archives-ouvertes.fr/hal-02391939 (2019).

[37] Y.-J. Lee and J. Xu, New formulations, positivity preserving discretizations and stability analysis for non-Newtonian flow models. *Comput. Methods Appl. Mech. Eng.* **195** (2006) 1180–1206.
[5] J.W. Barrett and S. Boyaval, Existence and approximation of a (regularized) Oldroyd-B model. *Math. Models Methods Appl. Sci.* **21** (2011) 1783-1837.

[6] J.W. Barrett and S. Boyaval, Finite element approximation of the FENE-P model. *IMA J. Numer. Anal.* **38** (2018) 1599-1660.

[7] J.W. Barrett and E. Süli, Existence of global weak solutions to some regularized kinetic models for dilute polymers. *Multiscale Model. Simul.* **6** (2007) 506-546.

[8] J.W. Barrett, C. Schwab and E. Süli, Existence of global weak solutions for some polymeric flow models. *Math. Models Methods Appl. Sci.* **15** (2005) 939-983.

[9] A. Bermudez and M.E. Vazquez, Upwind methods for hyperbolic conservation laws with source terms. *Comput. Fluids* **23** (1994) 1049-1071.

[10] R.B. Bird and J.M. Wiest, Constitutive equations for polymeric liquids. *Ann. Rev. Fluid Mech.* **27** (1995) 169-193.

[11] R. Botchorishvili, B. Perthame and A. Vasseur, Equilibrium schemes for scalar conservation laws with stiff sources. *Math. Comput.* **72** (2003) 131-157.

[12] F. Bouchut, H. Ounaissa and B. Perthame, Upwinding of the source term at interfaces for Euler equations with high friction. *Comput. Math. App.* **53** (2007) 361-375.

[13] S. Boyaval, Lid-driven-cavity simulations of Oldroyd-B models using free-energy-dissipative schemes, in Numerical Mathematics and Advanced Applications 2009. Springer (2010) 191-198.

[14] S. Boyaval, T. Lelièvre and C. Mangoubi, Free-energy-dissipative schemes for the Oldroyd-B model. *Math. Model. Numer. Anal.* **43** (2009) 523-561.

[15] S. Boyd and L. Vandenberghe, Convex Optimization. Cambridge University Press, Cambridge (2004).

[16] CALIF³S, A software components library for the computation of fluid flows. https://gforge.irsn.fr/gf/project/califs.

[17] G.-Q. Chen, C. D. Levermore and T.-P. Liu, Hyperbolic conservation laws with stiff relaxation terms and entropy. *Commun. Pure Appl. Math.* **47** (1994) 787-830.

[18] R. Comminal, J. Spangenberg and J.H. Hattel, Robust simulations of viscoelastic flows at high Weissenberg numbers with the streamfunction/log-conformation formulation. *J. Non-Newtonian Fluid Mech.* **223** (2015) 37-61.

[19] M. Crouzeix and P.A. Raviart, Conforming and nonconforming finite element methods for solving the stationary Stokes equations. *RAIRO Série Rouge* **7** (1973) 33-75.

[20] H. Damanik, J. Hron, A. Ouazzi and S. Turek, A monolithic FEM approach for the log-conformation reformulation (LCR) of viscoelastic flow problems. *J. Non-Newtonian Fluid Mech.* **165** (2010) 1105-1113.

[21] A. Duran, F. Marche, R. Turpault and C. Berthon, Asymptotic preserving scheme for the shallow water equations with source terms on unstructured meshes. *J. Comput. Phys.* **287** (2015) 184-206.

[22] R. Fattal and R. Kupferman, Constitutive laws for the matrix-logarithm of the conformation tensor. *J. Non-Newtonian Fluid Mech.* **123** (2004) 281-285.

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

[24] T. Gallouët, R. Herbin, J.-C. Latché and K. Mallem, Convergence of the Marker-And-Cell scheme for the incompressible Navier-Stokes equations on non-uniform grids. *Found. Comput. Math.* **18** (2018) 249-289.

[25] D. Grapsas, R. Herbin, W. Kheriji and J.-C. Latché, An unconditionally stable staggered pressure correction scheme for the compressible Navier-Stokes equations. *SMAI J. Comput. Math.* **2** (2016) 51-97.

[26] J.M. Greenberg and A.Y. Leroux, A well-balanced scheme for the numerical processing of source terms in hyperbolic equations. *SIAM J. Numer. Anal.* **33** (1996) 1-16.

[27] A.M. Grillet, B. Yang, B. Khomami and E.S. Shaqfeh, Modeling of viscoelastic lid driven cavity flow using finite element simulations. *J. Non-Newtonian Fluid Mech.* **88** (1999) 99-131.

[28] F.H. Harlow and A.A. Amsden, A numerical fluid dynamics calculation method for all flow speeds. *J. Comput. Phys.* **8** (1971) 197-213.

[29] F.H. Harlow and J.E. Welsh, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface. *Phys. Fluids* **8** (1965) 2182-2189.

[30] R. Herbin, W. Kheriji and J.-C. Latché, On some implicit and semi-implicit staggered schemes for the shallow water and Euler equations. *Math. Model. Numer. Anal.* **48** (2014) 1807-1857.

[31] D. Hu and T. Lelievre, New entropy estimates for the Oldroyd-B model and related models. *Commun. Math. Sci.* **5** (2007) 909-916.

[32] M.A. Hulsen, R. Fattal and R. Kupferman, Flow of viscoelastic fluids past a cylinder at high Weissenberg number: stabilized simulations using matrix logarithms. *J. Non-Newtonian Fluid Mech.* **127** (2005) 27-39.

[33] S. Jin, Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations. *SIAM J. Sci. Comput.* **21** (1999) 441-454.

[34] T. Katsaounis, B. Perthame and C. Simeoni, Upwinding sources at interfaces in conservation laws. *Appl. Math. Lett.* **17** (2004) 309-316.

[35] J.-C. Latché and K. Saleh, A convergent staggered scheme for the variable density incompressible Navier-Stokes equations. *Math. Comput.* **87** (2018) 581-632.

[36] J.-C. Latché, B. Piar and K. Saleh, A discrete kinetic energy preserving convection operator for variable density flows on locally refined staggered meshes. https://hal.archives-ouvertes.fr/hal-02391939 (2019).

[37] Y.-J. Lee and J. Xu, New formulations, positivity preserving discretizations and stability analysis for non-Newtonian flow models. *Comput. Methods Appl. Mech. Eng.* **195** (2006) 1180-1206.
[38] R.J. LeVeque and W. Jinghua, A linear hyperbolic system with stiff source terms, in *Nonlinear Hyperbolic Problems: Theoretical, Applied, and Computational Aspects: Proceedings of the Fourth International Conference on Hyperbolic Problems*, Taormina, Italy, April 3 to 8, 1992. Friedrich Vieweg & Sohn, Wiesbaden (1993) 401–408.

[39] E. Lorin and V. Seignole, Convection systems with stiff source terms. *Math. Models Methods Appl. Sci.* **13** (2003) 971–1018.

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

[41] M. Lukáčová-Medvid'ová, H. Notsu and B. She, Energy dissipative characteristic schemes for the diffusive Oldroyd-B viscoelastic fluid. *Int. J. Numer. Methods Fluids* **81** (2016) 523–557.

[42] L. Moreno, R. Codina, J. Baiges and E. Castillo, Logarithmic conformation reformulation in viscoelastic flow problems approximated by a VMS-type stabilized finite element formulation. *Comput. Methods Appl. Mech. Eng.* **354** (2019) 706–731.

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

[44] R. Rannacher and S. Turek, Simple nonconforming quadrilateral Stokes element. *Numer. Methods Part. Differ. Equ.* **8** (1992) 97–111.

[45] M. Renardy, Mathematical Analysis of Viscoelastic Flows. SIAM (2000).

[46] P. Saramito, On a modified non-singular log-conformation formulation for Johnson-Segalman viscoelastic fluids. *J. Non-Newtonian Fluid Mech.* **211** (2014) 16–30.

[47] M.F. Tomé, A. Castelo, A.M. Afonso, M.A. Alves and F. Pinho, Application of the log-conformation tensor to three-dimensional time-dependent free surface flows. *J. Non-Newtonian Fluid Mech.* **175** (2012) 44–54.

[48] K. Yapici, B. Karasoz and Y. Uludag, Finite volume simulation of viscoelastic laminar flow in a lid-driven cavity. *J. Non-Newtonian Fluid Mech.* **164** (2009) 51–65.

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

Please help to maintain this journal in open access!

This journal is currently published in open access under the **Subscribe to Open model (S2O)**. We are thankful to our subscribers and supporters for making it possible to publish this journal in open access in the current year, free of charge for authors and readers.

Check with your library that it subscribes to the journal, or consider making a personal donation to the S2O programme by contacting subscribers@edpsciences.org.

More information, including a list of supporters and financial transparency reports, is available at https://edpsciences.org/en/subscribe-to-open-s2o.
[38] R.J. LeVeque and W. Jinghua, A linear hyperbolic system with stiff source terms, in *Nonlinear Hyperbolic Problems: Theoretical, Applied, and Computational Aspects: Proceedings of the Fourth International Conference on Hyperbolic Problems*, Taormina, Italy, April 3 to 8, 1992. Friedrich Vieweg & Sohn, Wiesbaden (1993) 401–408.

[39] E. Lorin and V. Seignole, Convection systems with stiff source terms. *Math. Models Methods Appl. Sci.* **13** (2003) 971–1018.

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

[41] M. Lukáčová-Medvid'ová, H. Notsu and B. She, Energy dissipative characteristic schemes for the diffusive Oldroyd-B viscoelastic fluid. *Int. J. Numer. Methods Fluids* **81** (2016) 523–557.

[42] L. Moreno, R. Codina, J. Baiges and E. Castillo, Logarithmic conformation reformulation in viscoelastic flow problems approximated by a VMS-type stabilized finite element formulation. *Comput. Methods Appl. Mech. Eng.* **354** (2019) 706–731.

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

[44] R. Rannacher and S. Turek, Simple nonconforming quadrilateral Stokes element. *Numer. Methods Part. Differ. Equ.* **8** (1992) 97–111.

[45] M. Renardy, Mathematical Analysis of Viscoelastic Flows. SIAM (2000).

[46] P. Saramito, On a modified non-singular log-conformation formulation for Johnson-Segalman viscoelastic fluids. *J. Non-Newtonian Fluid Mech.* **211** (2014) 16–30.

[47] M.F. Tomé, A. Castelo, A.M. Afonso, M.A. Alves and F. Pinho, Application of the log-conformation tensor to three-dimensional time-dependent free surface flows. *J. Non-Newtonian Fluid Mech.* **175** (2012) 44–54.

[48] K. Yapici, B. Karasoz and Y. Uludag, Finite volume simulation of viscoelastic laminar flow in a lid-driven cavity. *J. Non-Newtonian Fluid Mech.* **164** (2009) 51–65.

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

Please help to maintain this journal in open access!

This journal is currently published in open access under the **Subscribe to Open model (S2O)**. We are thankful to our subscribers and supporters for making it possible to publish this journal in open access in the current year, free of charge for authors and readers.

Check with your library that it subscribes to the journal, or consider making a personal donation to the S2O programme by contacting subscribers@edpsciences.org.

More information, including a list of supporters and financial transparency reports, is available at https://edpsciences.org/en/subscribe-to-open-s2o.