Temporal Artificial Stress Diffusion for Numerical Simulations of Oldroyd-B Fluid Flow

Abstract

This paper presents a numerical evaluation of two different artificial stress diffusion techniques for the stabilization of viscoelastic Oldroyd-B fluid flows at high Weissenberg numbers. The standard artificial diffusion in the form of a Laplacian of the extra stress tensor is compared with a newly proposed approach using a discrete time derivative of the Laplacian of the extra stress tensor. Both methods are implemented in a finite element code and demonstrated in the solution of a viscoelastic fluid flow in a two-dimensional corrugated channel for a range of Weissenberg numbers. The numerical simulations have shown that this new temporal stress diffusion not only efficiently stabilizes numerical simulations, but also vanishes when the solution reaches a steady state. It is demonstrated that in contrast to the standard tensorial diffusion, the temporal artificial stress diffusion does not affect the final solution.

1. Introduction

The mathematical modeling and numerical simulation of non-Newtonian fluid flows became one of the fastest developing and most challenging problems of contemporary computational fluid dynamics. It has numerous physical and technical applications of high interest. One of the key applications motivating the below-described work comes from hemodynamics. Due to the complex microstructure of blood, it can be considered as a viscoelastic fluid, sharing some characteristic behaviors with polymeric liquids. The most distinct property of viscoelastic fluids is the ability to store and recover mechanical energy. Such fluids, including blood, can be mathematically described by the class of tensorial constitutive relations, from which the Oldroyd-B model (possibly generalized) is one of the simplest and most often used [1,2].

1.1. Modeling Challenges

Due to the complex behavior of viscoelastic fluids and their rheological description, mathematical modeling is a rather challenging problem. The classical set of balance equations (usually represented by Navier–Stokes equations) has to be supplemented by suitable tensorial constitutive relation between the stress and the rate of strain or rate of deformation.

One of many typical mathematical problems in the modeling of viscoelastic fluids is related to the existence of multiple characteristic time scales of the fluid flow. These so-called relaxation/retardation times in some sense allow these fluids to remember (or to feel) their deformation and stress history. Therefore, the model is more sensitive to the prescribed initial and boundary data, and extra attention should be paid to their choice and numerical treatment.

From the computational point of view, apart from the additional equations to discretize more complex and sensitive computational setups, the biggest challenge arises with large relaxation times. This issue is manifested by a breakdown of stability numerical solutions for larger relaxation time values. In the literature, this problem is extensively reported and studied, referred to as the high Weissenberg number problem, where the Weissenberg number is an essential non-dimensional parameter (proportional to the relaxation time) of viscoelastic fluid flows.

To deal with this problem, numerous approaches have been proposed and tested by various authors. One of the commonly used techniques to overcome or at least minimize the numerical instabilities arising at high Weissenberg numbers is the addition of extra artificial stabilization to the numerical solver. There is a wide body of specialized literature describing and dealing with this issue [3,4,5,6,7,8,9,10,11,12,13,14].

One of the simplest and most frequently used stabilization methods is based on the addition of an artificial stress diffusion term into the constitutive equations for a viscoelastic stress tensor. Such an additional tensorial term can be seen as an artificial diffusion applied separately to individual stress components. Although this numerical diffusion is usually used as a purely artificial term, there are also physical arguments, based on the microstructural concept of fluids rheology, for the addition of such a physical term into the constitutive model [15,16,17]. Independently of interpretation, the additional term modifies the rheological constitutive model and evidently affects (alters) the solution of the considered problem. This is why attention should be paid to keep the additionally modified model consistent with (or at least close enough to) the original problem [18,19].

1.2. Goal and Structure of the Paper

The goal of this paper is to present and test a new temporal stress diffusion technique that can avoid the unwanted artificial artifacts of the standard numerical stress stabilization.

In place of the widely used artificial stress diffusion based on the Laplacian of extra stress $\alpha \Delta \mathit{\tau }$, a newly proposed term corresponding to the Laplacian of (pseudo) time derivative of extra stress $\alpha \Delta {\mathit{\tau }}_t$ is used and tested here.

In this paper, we first show in detail how the standard stress diffusion stabilization affects the solution for a range of Weissenberg numbers. Then, we introduce and test a novel temporal artificial stress diffusion that preserves the stabilization capabilities of the former standard technique, while avoiding the unwanted and non-physical alteration of the final solution. Both these properties of the new technique are demonstrated in a number of simulations. The main contribution and novelty of our work can thus be seen in the introduction and testing of this new temporal stress diffusion technique for the stabilization of viscoelastic fluid flow simulations.

The structure of the paper is as follows. Starting from the introductory Section 1, the mathematical model is summarized in Section 2, and the numerical method is described in Section 3. The two (standard and new) tensorial diffusion stabilization methods are described in Section 3.2. Numerical simulations are summarized in Section 4, where, first, the standard stabilization is presented, followed by a similar set of simulation results obtained using the newly proposed temporal stabilization. All the results are discussed and the methods are compared. A final summary of findings and some suggested extensions and outlooks are presented in Section 5.

This paper builds on our first short study published in conference proceedings [20]. Said work briefly presented the temporal stress diffusion principle and showed that it may have a positive impact on the stability of the numerical algorithm. The few initial results presented in this conference contribution were promising and showed a mild dependence of the stabilizing effect on the Reynolds number. This initial success encouraged a much larger and deeper study presented in full in this paper. Here, we have tried to present an extended study of the new temporal stress diffusion stabilization. We kept the geometry of the test case, but this time we worked with a fixed Reynolds number (typical for blood flows) and compared the effects of two methods, the standard and the new temporal stress diffusion, on the stability of the numerical algorithm for a range of Weissenberg numbers. The results presented herein are based on a comprehensive set of new numerical simulations performed to verify and support the conclusions of this paper.

2. Mathematical Model

Governing Equations

The following non-dimensional governing equations, defined in a bounded domain $\mathsf{\Omega }\subset {\mathbb{R}}^d(d=2,3)$ over the time interval $\left[0,T_f\right]$, model the unsteady, incompressible, isothermal viscoelastic flow of homogeneous Oldroyd-B fluid considered in this work:

$$\left\{\begin{array}{c}\mathcal{R}e\left({\displaystyle \frac{\partial \mathit{u}}{\partial t}}+\mathit{u}·\nabla \mathit{u}\right)+\nabla p=2(1-\eta )\nabla ·\mathsf{D}+\nabla ·\mathit{\tau }\hfill \\ \nabla ·\mathit{u}=0\hfill \\ \mathit{\tau }+\mathcal{W}e\left({\displaystyle \frac{\partial \mathit{\tau }}{\partial t}}+\mathit{u}·\nabla \mathit{\tau }-\nabla {\mathit{u}}^T·\mathit{\tau }-\mathit{\tau }·\nabla \mathit{u}\right)=2\eta \mathsf{D}\hfill \end{array}\right.$$

(1)

Here, $\mathit{\tau }$ is the viscoelastic (or polymeric) part of the complete stress tensor $\mathsf{T}$ = $2{\nu }_s\mathsf{D}+\mathit{\tau }$, where the symmetric part of the velocity gradient (rate of strain) tensor is denoted by $\mathsf{D}$ = $\frac{1}{2}\left(\nabla \mathit{u}+\nabla {\mathit{u}}^T\right)$. The viscosity ${\nu }_s$ usually represents the physical viscosity of the solvent (for a polymeric solution). Here, it is used to denote a part of the total (apparent) stress viscosity depending on the chosen stress model. In the case described later on, $\eta \in [0,1]$ is defined in such way that ${\nu }_s+\eta =1$ is the total adimensional kinematic viscosity. This means that when $\eta =0$, the governing system (1) is reduced to the Navier–Stokes equations for Newtonian fluid. If $\eta =1$, the system (1) reduces to the upper-convected Maxwell fluid model.

In contrast to the standard incompressible Newtonian fluid, the flow of a viscoelastic fluid is characterized by two essential non-dimensional parameters. The first one (known from Newtonian fluid mechanics) describes the relative effects of inertial and viscous forces and is given by their mutual rates, leading to the well known Reynolds number $\mathcal{R}e=\frac{\mathit{U}L}{\nu }$. The additional non-Newtonian parameter is the Weissenberg number$\mathcal{W}e=\frac{{\lambda }_1\mathit{U}}{L}$, characterizing the elasticity of the fluid. One of the possible interpretations is to consider this non-dimensional parameter as the rate of two distinct time scales, one of the memory of the fluid (described by the relaxation time ${\lambda }_1$) and the second of the convection/advection time scale $L/\mathit{U}$, indicating the time needed by the fluid moving at the characteristic velocity U to pass the characteristic distance L.

3. Numerical Solution

The details of the discretization of the model (1), via the characteristic finite element method, can be found in our previous work [20]. Here, we just recall the semi-discrete variational form of the system, where the time advancing part of the algorithm, needed to explain the stabilization technique, can be seen.

3.1. Semi-Discrete Variational Form

The finite set of time instants $t_n=n\Delta t=n\frac{T_{\mathit{f}}}{N},$ with n = $0,$ $\dots ,N$ is defined over the given interval $[0,T_{\mathit{f}}]$.

Both the momentum and constitutive equations are discretized with respect to time using the implicit Euler scheme and characteristic Galerkin method for the convective terms.

For the equations of motion (momentum and continuity balance), it leads to

$$\left\{\begin{array}{cc}{\displaystyle {\int }_{\mathsf{\Omega }}2(1-\eta ){\mathsf{D}}^n:\nabla \mathit{v}+\mathcal{R}e{\int }_{\mathsf{\Omega }}{\displaystyle \frac{{\mathit{u}}^n-{\mathit{u}}_{★}^{n-1}}{\Delta t}}·\mathit{v}-{\int }_{\mathsf{\Omega }}{p}^n\nabla ·\mathit{v}={\int }_{\mathsf{\Omega }}\left(\nabla ·\mathit{\tau }\right)\mathit{v},}\hfill & \forall \mathit{v}\in {\mathbf{H}}_0^1\left(\mathsf{\Omega }\right)\hfill \\ {\displaystyle {\int }_{\mathsf{\Omega }}\left(\nabla ·{\mathit{u}}^n\right)q=0,}\hfill & \forall q\in L_0^2\left(\mathsf{\Omega }\right)\hfill \end{array}\right.$$

(2)

while for the discretization of tensorial constitutive equation it gives

$${\displaystyle {\int }_{\mathsf{\Omega }}\left({\mathit{\tau }}^n+\mathcal{W}e{\displaystyle \frac{{\mathit{\tau }}^n-{\mathit{\tau }}_{★}^{n-1}}{\Delta t}}\right):\mathit{S}={\int }_{\mathsf{\Omega }}\left[2\eta \mathsf{D}+\mathcal{W}e\left(\nabla {\mathit{u}}^T·{\mathit{\tau }}^n+{\mathit{\tau }}^n·\nabla \mathit{u}\right)\right]:\mathit{S},\forall \mathit{S}\in \mathcal{S}}$$

(3)

More details on the implementation, convergence, stability and other properties of this specific semi-discretization can be found in our previous paper [20] or in references [21,22].

3.2. Stabilization–Stress Diffusion

The numerical stabilization for the above-described algorithm under a high Weissenberg number regime was carried out using an additional tensorial diffusive term $\mathit{E}$ introduced to the elastic stress constitutive equation

$$\mathit{\tau }+\mathcal{W}e\left({\displaystyle \frac{\partial \mathit{\tau }}{\partial t}}+\mathit{u}·\nabla \mathit{\tau }-\nabla {\mathit{u}}^T·\mathit{\tau }-\mathit{\tau }·\nabla \mathit{u}\right)=2\eta \mathsf{D}+\mathit{E}.$$

In this work, we will use and compare two different versions of the additional term $\mathit{E}$.

• The standard diffusive term proportional to the Laplacian of elastic stress:

  $$\mathit{E}=\alpha ·\Delta \mathit{\tau }\approx \alpha ·\Delta {\mathit{\tau }}^n$$

  (4)
  This standard diffusion has frequently been used by several researchers. In general, this extra term does not vanish when the solution reaches the steady state, so $\mathit{E}=\alpha ·\Delta {\mathit{\tau }}^n\nrightarrow 0$ when ${\mathit{\tau }}^n\stackrel{n\to \infty }{\to }\mathit{\tau }$. This means that the original problem is permanently modified by the added diffusive term and thus results may be quite sensitive to the values of the parameter $\alpha$.
• The new temporal diffusive term is proportional to the Laplacian of the (pseudo) time derivative of the elastic stress:

  $$\mathit{E}=\alpha ·\Delta {\mathit{\tau }}_t\approx \alpha ·\Delta \left({\mathit{\tau }}^n-{\mathit{\tau }}^{n-1}\right)$$

  (5)
  Here, the temporal index n corresponds to pseudo-time, used in the time-marching iterative procedure. In this case, the extra term $\mathit{E}$ will vanish when the numerical solution converges to the steady state ${\mathit{\tau }}^n\to \mathit{\tau }$, i.e., $\left({\mathit{\tau }}^n-{\mathit{\tau }}^{n-1}\right)\to \mathbf{0}$. Due to this property, the added diffusivity should only act during the initial, transitional stage of (pseudo) time stepping. Thus, the solution of the original model should be recovered (the stress diffusion will vanish) and the final results should not be sensitive to the choice of the parameter $\alpha$.

4. Numerical Simulations

As a first step in the analysis of the efficiency of the new temporal diffusive term, we consider a steady flow test case case. For simplicity, keeping in mind the future unsteady applications, the steady problem solution is searched by the time-marching technique, i.e., by solving an unsteady problem, where the steady state is recovered for $t\to \infty$, subject to stationary boundary conditions.

The in-house written code (using FreeFem++ toolboxes [23]) was first validated by comparing the numerical solution of the Poiseuille flow in the straight pipe (2D channel) with the corresponding analytical solution for the Oldroyd-B model, which can be found, e.g., in [1] (this validation is not presented here). Within this paper, a more complex tube-like geometry modification is adopted. A multiply narrowed (corrugated) channel was considered to demonstrate the effects of numerical stabilization for non-trivial flows at higher Weissenberg numbers. One of the features that led us to the choice of this non-standard test geometry (already used in [6]) is that, even for stationary boundary conditions, the stress acting on a flowing fluid parcel has a quasi periodical nature. The repeatedly contracting and expanding channel cross section successively acts by an elevated local stress to the fluid passing along the channel. Such successive loading–unloading cycles are of interest in studying viscoelastic fluids flows, where the time history of stress plays an important role.

The numerical simulations presented hereafter were performed both with and without the additional numerical stabilization. Two different versions of the stabilization term $\mathit{E}$ were considered and several values of stabilization parameters $\alpha$ were used.

The (dimensionless) elastic viscosity $\eta =0.1$ was used in all tests. The Reynolds number was kept fixed at the value $\mathcal{R}e=1000$. The Weissenberg number was increased incrementally, while taking the previous converged solution as initial condition for higher $\mathcal{W}e$ simulation. This approach is often referred to as the continuation method.

The standard flow field variables (fluid velocity, pressure and tensor fields) were evaluated in all numerical simulations. In addition, the elastic stress tension on the wall was computed as

$${\mathit{\tau }}_w=-(\mathit{\tau }·\widehat{\mathit{n}})·\widehat{\mathit{t}}{|}_w$$

(6)

• Boundary conditions

The straight (non-corrugated) inlet and outlet parts of the channel allow us to consider the fully developed flow in the vicinity of the respective boundaries and thus it can be used to set the boundary conditions. The Dirichlet boundary condition on the inlet was determined by the analytical Poiseuille solution for the Oldroyd-B fluid model. The standard no-slip (i.e., homogeneous Dirichlet) conditions were imposed for the fluid velocity.

• Computational domain

The 2D axi-symmetrical sinusoidally corrugated (narrowed) pipe was assumed; see Figure 1. The contracting/expanding parts have lengths of $L_{con}$ and $L_{exp}$. Between these parts, several identical segments characterized by diameters $D_{min}$ resp. $D_{max}$ and the length $L_{seg}$ were inserted.

The near regularity shown in Figure 2 was generated by FreeFem++, using the Delaunay–Voronoi algorithm. The size of mesh cells was set by dividing the boundary into 10 equal segments per unit of length. No special treatment was used to symmetrize the grid with respect to the x axis.

Table 1. Characteristic parameters of the mesh.

Elements	Nodes ${\mathit{P}}_2$	Nodes ${\mathit{P}}_1$	${\mathit{h}}_{\mathit{max}}$	${\mathit{h}}_{\mathit{min}}$
3072	6485	1707	0.0844012	0.0395281

provides some data on the mesh used.

4.1. Numerical Results

The viscoelastic fluid flow through the 2D corrugated channel was solved for a range of Weissenberg numbers at a single fixed Reynolds number. The aim of the performed numerical simulations was to find the limits of applicability, in terms of maximum attainable Weissenberg number $\mathcal{W}e$, of the original code (without artificial diffusion) and its updated version, using both versions of the tensorial artificial diffusion terms (4) and (5).

For the above-described test case, using only the reference viscosity scheme (without added tensorial diffusion), the maximum Weissenberg number $\mathcal{W}e\approx 0.4$ was reached. Using the same algorithm with additional classical stress diffusion (4), with constant or variable diffusion parameter $\alpha$, the maximum attained Weissenberg number can be significantly higher. However, depending on the amount of artificial stress diffusion (value of $\alpha$), the obtained numerical solution may be quite far from the numerical solution of the original Oldroyd-B fluid flow problem.

4.1.1. Numerical Solutions for the Standard Artificial Diffusion

Several test were performed with the standard stress diffusion $\mathit{E}=\alpha ·\Delta {\mathit{\tau }}^n$ considering different (three constant and one variable) values of parameter $\alpha \in \{{10}^{-4},{10}^{-3},{10}^{-2},\alpha \left(\mathcal{W}e\right)\}$.

The idea of making the numerical diffusion coefficient $\alpha$ variable, depending on Weissenberg number, is motivated by the fact that for low $\mathcal{W}e$, no stabilization is required, but for higher and increasing values of $\mathcal{W}e$, the numerical stabilization becomes necessary, so the coefficients of the numerical diffusion have to increase. On the other hand, the added diffusion modifies the final solution and thus it should be kept as small as possible. Thus, the aim was to make the adjustment of suitable numerical diffusion coefficients automatic, always staying close to the minimum necessary value needed to stabilize the numerical solution, but avoiding an over-stabilized i.e., over-smoothed solution.

Based on the above-mentioned goal, it is clear that to set the numerical diffusion coefficient one needs a function (depending on Weissenberg number) that will be monotonically increasing with $\mathcal{W}e$, having very small values up to certain threshold, after which it will grow up smoothly up to certain asymptotic value for very high Weissenberg numbers. The existence of the threshold, below which (almost) no numerical diffusion is needed, was observed in our preliminary simulations (and the corresponding critical $\mathcal{W}e$ was known). The need for an upper asymptotic value comes from the fact that the added diffusion significantly modifies the model and its results, and thus for very high added diffusion, the obtained results are no closer to the solution of the original non-diffusive problem.

These requirements lead to the choice of the function arctan, which is smooth, increasing with values asymptotically bounded from below and above, with a clearly defined rapid increase region. This function was used and properly shifted and scaled. The proportionality coefficient $\Delta t{h}^2$ was added to keep the discrete model consistent with the continuous one, guaranteeing that the artificial diffusion term will vanish as the temporal step $\Delta t$ and spatial step h tend to zero.

The variable parameter $\alpha \left(\mathcal{W}e\right)$ depends on the Weissenberg number as $\alpha \left(\mathcal{W}e\right)\propto h^2\Delta t·\mathrm{atan}\left(\epsilon \mathcal{W}e\right)$, with asymptotic values ${\alpha }_0=0.0$ and ${\alpha }_{\infty }=0.04$ (see Figure 3).

Figure 4 shows the graphs of the norm of the stabilization term $\parallel \mathit{E}\parallel =\parallel \alpha ·\Delta {\mathit{\tau }}^n\parallel$ with respect to the number of iterations for three different values of the stabilization parameter $\alpha$. From this figure several important observations can be made.

• When increasing the value of $\alpha$, the value of $\parallel \mathit{E}\parallel =\parallel \alpha ·\Delta {\mathit{\tau }}^n\parallel$ increases as well, but not directly proportionally to $\alpha$. The value of $\alpha$ increases 100 times between the left and right plots, but the norm $\parallel \mathit{E}\parallel =\parallel \alpha ·\Delta {\mathit{\tau }}^n\parallel$ only increases about 10 times. This can be attributed to the smoothing effect of the diffusive term $\mathit{E}$ leading to the reduction in the norm of the Laplacian of the elastic stress tensor $\parallel \Delta {\mathit{\tau }}^n\parallel$.
• The (norm of the) added extra diffusive term $\mathit{E}$ never vanishes, so the converged solution corresponds to another (modified) problem other than the original one without the diffusive term.
• With higher values of the stabilization parameter $\alpha$, a solution can be obtained for a higher Weissenberg number (graph lines with different colors correspond to different $\mathcal{W}e$). However, at the same time the norm of the additional term $\parallel \mathit{E}\parallel$ grows with $\alpha$, so the obtained solution is further from the sought solution of the original, non-diffusive Oldroyd-B fluid flow model.
• For higher values of $\alpha$ and a higher-value Weissenberg number, many more iterations are needed to achieve the steady, converged solution.

As mentioned before, when trying to improve the behavior of the standard diffusive term, a variable coefficient $\alpha \left(\mathcal{W}e\right)$ was tested. The dependence of $\alpha$ on the Weissenberg number was set in such a way that for small Weissenberg numbers, the $\alpha$ was also small and only increased for elevated Weissenberg numbers, when the instability occurs. The results for this variable $\alpha \equiv \alpha \left(\mathcal{W}e\right)$ can be seen in Figure 5, leading to observations similar to those for the constant $\alpha$. In this variable $\alpha$ case, it is possible to obtain stable steady convergent solutions for $\mathcal{W}e=3.76$.

The summary of the highest attainable Weissenberg numbers for given test case and different settings of $\alpha$ can be seen in

Table 2. Maximum value of the Weissenberg number for which the steady numerical solution is obtained, depending on $\alpha$.

$\mathit{\alpha }$	0	${10}^{-4}$	${10}^{-3}$	${10}^{-2}$	$\mathit{\alpha }\left(\mathcal{W}\mathit{e}\right)$
${\mathcal{W}e}_{max}$	0.44	0.63	1	2	3.76

.

The added standard stabilization $\mathit{E}=\alpha ·\Delta \mathit{\tau }$ really leads to an increase in the maximum attainable Weissenberg number. However, it should be stressed again that the added tensorial diffusion modifies the original problem and thus the numerically found solution does not correspond (especially for higher $\alpha$) to the original (non-diffusive) Oldroyd-B model.

The smoothing effect introduced by the standard diffusive term $\mathit{E}$ can be assessed from the Figure 6 showing the norm of the second derivative (Laplacian) of the elastic tensor for the numerical solutions obtained without stabilization ($\alpha =0$) and with stabilization ($\alpha >0$). The values of the Laplacian norm decrease with the increase in the diffusion parameter $\alpha$. This demonstrates the fact that each such numerical solution corresponds to different diffusive problem. Even for very small values of $\alpha$, the standard stabilization significantly modifies the solution in comparison with the classical Oldroyd-B problem.

The nature of the solution and its modification by the added standard stress diffusion can be seen from the following Figure 7, Figure 8 and Figure 9, showing the contour maps of the components of the extra stress tensor $\mathit{\tau }$ obtained for the same case (moderate Weissenberg number) for different values of $\alpha$.

The viscoelastic stress tensor $\mathit{\tau }$ has the form

$$\mathit{\tau }=\left(\begin{array}{cc}{\tau }_{11}& {\tau }_{12}\\ {\tau }_{12}& {\tau }_{22}\end{array}\right),$$

(7)

so in the chosen coordinate system the ${\tau }_{11}$ corresponds to axial, ${\tau }_{12}$ to shear and ${\tau }_{22}$ normal component of the tensor.

From the shown tensor fields is obvious that even the smallest values of $\alpha$ lead to visible (and non-negligible) differences between the solutions of the original (non-diffusive) and the modified Oldroyd-B model with artificial stress diffusion.

For example, the axial stress component ${\tau }_{11}$ shown in the Figure 7 exhibits visible peaks (local maxima) in the near wall layer in the contractions in the non-diffusive case (for $\alpha =0$), but these peaks are largely suppressed and smoothed out in all the cases where the artificial diffusive stress term was used. From the sequence of results for $\alpha ={10}^{-4},{10}^{-3},{10}^{-2}$ is clear that higher coefficient $\alpha$ leads to more pronounced smoothing effect. Even the smallest value $\alpha ={10}^{-4}$ leads to a tensor field where the stress maxima are visibly cut off, although the overall nature of the field remains similar. This behavior and dependence on the diffusion coefficient $\alpha$ is also confirmed for the shear stress component ${\tau }_{12}$ shown in Figure 8 and the normal stress component ${\tau }_{22}$ shown in Figure 9.

It is worth noting that, judging from Figure 7, Figure 8 and Figure 9, the choice of diffusion coefficient $\alpha \left(\mathcal{W}e\right)$ depending on the Weissenberg number leads to highly smoothed results that are quite far from those obtained without an artificial diffusive term.

The smoothing effect of the standard stress diffusion can also be assessed considering the plots of extremal values of individual stress components shown in Figure 10. The curves marked with $\alpha =0$ correspond to the sought solution of the (non-diffusive) Oldroyd-B model, while the others represent the results obtained using a different amount of stress diffusion, of the standard kind. Evidently, when using the standard stress diffusion, the solution can be found for higher Weissenberg numbers, but the values of the stress extrema are very different (smaller) than in the case solved without the stress diffusion. It should be noted that the variable $\alpha \left(\mathcal{W}e\right)$ slightly improves the results for smaller Weissenberg numbers, but for higher $\mathcal{W}e$, it set quite high values of $\alpha$ (see the numbers along the corresponding lines in the plot), which leads to over smoothed solution.

One of the most important global flow parameters, with high practical importance, is the pressure drop. It is defined as the difference in total pressure between the two ends of the pipe (channel). It occurs when frictional forces, caused by the resistance to flow, act on a fluid as it flows through the tube. The main parameters affecting the resistance of fluid to flow are the fluid velocity and viscosity (constant in this case). Pressure drop increases proportionally to the frictional shear forces within the tube. It can be partially assessed by comparing the fields of stress components in the Figure 10 (middle) and Figure 11.

Higher local flow velocities result in a larger pressure drop, while a low velocity will result in a lower or no pressure drop. From Figure 11, it is possible to confirm that the pressure drop decreases with the added stress diffusion. This behavior is qualitatively the same for all values of parameter $\alpha$.

Due to extra smoothing, introduced by the stress diffusion, the elastic stress tension on the wall ${\mathit{\tau }}_w$ (defined by (6)) decreases when the diffusion parameter $\alpha$ grows for a given Weissenberg number, as seen in Figure 12. They also show that the increase in $\mathcal{W}e$ leads to appropriate increase in the elastic stress and thus also its components on the wall. The stress diffusion term tends to symmetrize the graphs of the tension, while the elastic stress tension on the wall is naturally non-symmetric for the classical (non-diffusive) Oldroyd-B problem.

The results obtained using the standard stress diffusion $\mathit{E}=\alpha ·\Delta {\mathit{\tau }}^n$ have clearly shown that the method becomes numerically stable also for higher Weissenberg numbers. However the stability comes at the expense of higher necessary diffusion parameters $\alpha$, leading to a solution that might be quite different from the one corresponding to the original (non-diffusive) Oldroyd-B fluid flow model. The results presented in this section should be directly compared with those in the next section, showing the new temporal stress diffusion effects.

4.1.2. Numerical Solutions for the New Temporal Stress Diffusion

The same set of test simulations as for the standard stress diffusion were also performed for the new, temporal stabilization of the form $\mathit{E}=\alpha ·\Delta \left({\mathit{\tau }}^n-{\mathit{\tau }}^{n-1}\right)$. Additionally, for this modified version of stress diffusion several values of parameter $\alpha$ were considered, $\alpha \in \{{10}^{-4},{10}^{-3},{10}^{-2},{10}^{-1},\alpha \left(\mathcal{W}e\right)\}$, where the variable function $\alpha \left(\mathcal{W}e\right)$ was described in the Section 4.1.1. A summary is shown in the Figure 13.

The results obtained using the new temporal stress diffusion stabilization are quite different to those obtained using the standard stress diffusion. In contrast to all the previous results, the simulations shown here for the same Weissenberg number, the norm of the stabilization term $\alpha ·\Delta \left({\mathit{\tau }}^n-{\mathit{\tau }}^{n-1}\right)$, is almost insensitive to the diffusive parameter $\alpha$. Although in some cases more iterations are needed to reach a steady solution, the norm of the new temporal artificial stress diffusion is much smaller, and vanishing in time, when compared with the norm for the standard stabilization term (shown in Figure 4, Figure 5 and Figure 14).

This insensitivy to the parameter $\alpha$ is also reflected in the maximum attainable value of the Weissenberg number for which it is possible to obtain a steady numerical solutions. It is almost the same, around 0.68, independently of $\alpha$ (see

Table 3. Maximum value of Weissenberg number for which the steady numerical solution is obtained, depending on $\alpha$.

$\mathit{\alpha }$	0	${10}^{-4}$	${10}^{-3}$	${10}^{-2}$	${10}^{-1}$	$\mathit{\alpha }\left(\mathcal{W}\mathit{e}\right)$
${\mathcal{W}e}_{max}$	0.44	0.56	0.62	0.69	0.68	0.68

).

The most important observation, however, is that the numerical solutions obtained using the new temporal stress diffusion stabilization (5) are not affected by the stabilization term and all such solutions, independently of the parameter $\alpha$ used, correspond to the original classical Oldroyd-B model (without artificial diffusion).

The fields of individual stress components, the axial ${\tau }_{11}$, shear ${\tau }_{12}$ and normal ${\tau }_{22}$ are plotted in the Figure 15, Figure 16 and Figure 17, respectively. The results obtained without the artificial stress diffusion, i.e., with $\alpha =0$ are always shown on top of the sequence of contours for reference. From the visual comparison of tensor fields obtained for parameter $\alpha$ set to ${10}^{-4},{10}^{-3},{10}^{-2},{10}^{-1},$ is clear that the results are identical, showing no visible differences for various values of $\alpha$, including the reference solution with $\alpha =0$. This just demonstrates how robust the new temporal stabilization technique is, providing the required stabilization effect without unwanted artifacts. The temporal stress diffusion is thus very tolerant to the choice of the parameter $\alpha$ in (5), which can now only affect the convergence speed (number of iterations), but not the final result.

Such behavior was expected and desired, because the steady state was searched (and found), and the time derivatives of all quantities vanished and thus also the artificial diffusion term, proportional to time-derivative of the stress tensor, finally disappeared, without having a chance to affect the final solution. Considering this behavior, the automatic setting of $\alpha =\alpha \left(\mathcal{W}e\right)$ seems to be an easy and safe choice, because it eliminates the need for a user manual adjustment of the diffusion coefficient. This is confirmed in our simulations, where the first and last contour fields (corresponding to $\alpha =0$ and $\alpha =\alpha \left(\mathcal{W}e\right)$) in Figure 15, Figure 16 and Figure 17 are identical.

As for the standard stabilization, here some additional important information can be found in the plots of the maxima of the stress components with respect to Weissenberg number shown in Figure 18.

Obviously, both the non-stabilized ($\alpha =0$) and stabilized ($\alpha >0$) algorithms provide identical solutions up to the Weissenberg number $\mathcal{W}e=0.4$, above which the non-diffusive algorithm fails to converge, while the one with the new artificial temporal stress diffusion defined by (5) keeps providing stable solution up to $\mathcal{W}e$ ≈ $0.7$, as we can see in Figure 18.

Some limitations can be seen in adjusting the appropriate value of $\alpha$ for higher Weissenberg numbers. Taking into account the profile obtained for the case $\alpha ={10}^{-1}$, it seems to be possible to improve the adjustment of the variable function $\alpha \left(\mathcal{W}e\right)$, to reach stable solution for even higher Weissenberg numbers.

Using the new temporal artificial stress diffusion, the obtained pressure drop is completely comparable (see Figure 19 and Figure 20) with the results obtained for the original (non-diffusive) Oldroyd-B fluid model. Not only that there is a visible increase in the maximum attainable Weissenberg number when the artificial diffusion is employed, but as before, the solutions with and without stabilization are identical in the range where both algorithms converge. This agreement between the two algorithms is manifested in the Figure 19 as well as Figure 20, where both solutions are plotted.

In fact, the pressure drop does not changes with the diffusion parameter $\alpha$ (see Figure 20).

Consequently, the elastic stress tension on the wall ${\mathit{\tau }}_w$ does not changes with the parameter $\alpha$ and increases with Weissenberg number, in the same way as the solution of the classical non-diffusive Oldroyd-B model. Obviously, the temporal artificial stress diffusion term does not introduce any changes to the tension on the wall. The profile of the tension along the wall is plotted in Figure 21 for Weissenberg number $\mathcal{W}e=0.1$, showing the reference solution without any stabilization term compared with the solution obtained using the new temporal artificial stress diffusion. Both solutions apparently coincide, proving that the artificial extra stabilization term really vanished in the steady state.

For comparison, the same profile of the tension along the wall is shown in Figure 22 but for higher Weissenberg number $\mathcal{W}e=0.4$. In this case the profile shows higher, more sharp peaks (compared to $\mathcal{W}e=0.1$), with large negative values. This might be one of the reasons why for higher Weissenberg numbers the numerical simulations tend to become unstable, requiring the additional stabilization to converge.

The dramatic change in the stress tension on the wall between $\mathcal{W}e=0.1$ and $\mathcal{W}e=0.4$ is even more apparent in Figure 23, showing the corresponding profiles for different settings of the artificial diffusion parameter $\alpha$. This comparison again proves that the final result (obtained using new temporal stress diffusion) is independent of the choice of $\alpha$, which is in contrast to the behavior of the standard artificial diffusion, for which the results heavily depend on the choice of the diffusion parameter, as shown in Figure 12.

The main advantage and reason to use the new temporal artificial tensorial diffusion is that it makes possible to obtain the solution for higher values of the Weissenberg number, where the classical algorithm (without stabilization) easily becomes unstable, unable to converge to the steady state. For the artificial stress diffusion, proportional to Laplacian of (pseudo) time derivative of extra stress, not only that there is a visible increase in the maximum attainable Weissenberg number, when the artificial diffusion is employed, but as it was shown before, the solutions are identical in the range where both algorithms converge. This agreement between the two algorithms is apparent in Figure 18, Figure 20 and Figure 23, where both (stabilized and non-stabilized) solutions are plotted as one figure.

5. Conclusions and Remarks

The results presented in this paper significantly extend and deepen our knowledge concerning the new temporal stabilization method. The few initial simulations are already shown and discussed in [20], so the below-provided summary starts from where the conclusions in [20] ended.

• The numerical stability and robustness of the code were significantly improved by using an artificial diffusive stress term in the form (4), (5), respectively. This extra diffusive term allows one to obtain the solution of the given problem for higher values of $\mathcal{W}e$.
• It was shown that the standard tensorial diffusion (4) stabilizes the numerical method; however, it always affects the solution. So, strictly speaking, the obtained solution does not corresponds to the original Oldroyd-B fluid flow model. This means that this standard stabilization should be used with extreme caution and when it is used, the diffusion coefficient $\alpha$ should be as small as possible. It should be kept in mind that the same diffusive effect is also introduced by numerical diffusion contained in some low-order or upwinded schemes. Such methods should also be avoided or used with extreme caution.
• The newly proposed and tested temporal tensorial diffusion (5) has proved to be a simple and efficient method in solving high Weissenberg problems. For the given case, the maximum attainable Weissenberg number was increased by about 70% with respect to the original non-stabilized method. A further extension of this range seems to be possible by optimizing the choice of the (constant or variable) parameter $\alpha$.
• The main advantage of the new temporal diffusion is that it naturally vanishes when the solution converges towards the steady state. Thus, the final results are independent of the choice of the diffusion coefficient $\alpha$ and apparently the converged solution always corresponds to the solution of the original non-diffusive problem. In some sense, this temporal artificial diffusion resembles the concept of vanishing viscosity solutions to the (inviscid) Euler equations of gas dynamics.
• Another possible interpretation, or rather an analogy, of the proposed temporal vanishing stress diffusion is the residual smoothing technique. In case we interpret the difference of a certain quantity (e.g., stress tensor component in our case) between two consecutive pseudo time levels as a steady residual, then applying a Laplacian to it will have a similar effect like in the residual smoothing method that was largely popular for improving the numerical stability and convergence in CFD.
• Our future work will extend the verification of this type of temporal artificial diffusive term when using other numerical methods (and mathematical models), like the finite volumes and finite differences.
• A similar vanishing artificial diffusion effect can probably also be achieved by directly making the diffusion coefficient $\alpha$ vanish in time or taking it dependently on the time derivative of stress (or another quantity). Some initial tests for such an approach were shown recently in our proceedings paper [24].
• The question of the use of this new stabilization technique for unsteady problems remains open. It seems to be possible to consider it at least in the internal subiterations during the physical time stepping.
• The well-posedness of the complete continuous unsteady problem including the added Laplacian of the time derivative of the stress tensor is another very interesting and important open consideration for future investigations. It is far beyond the scope of our investigation, but there is a hope that such proof will be provided by someone in the future, considering that the added stabilization term is linear and particularly simple. Some hints can probably be found in older works related to certain pressure stabilization and projection methods for incompressible Navier–Stokes equations or in the analysis of some residual smoothing stabilization methods.

In conclusion, our simulations confirm that the new temporal artificial stress diffusion allows us to extend the range of usability of the numerical solver for higher Weissenberg numbers. For the considered numerical method and test case, the increase in the critical Weissenberg number was about 40%. From our experience, the variant of the stabilization algorithm with an automatic setting of $\alpha =\alpha \left(\mathcal{W}e\right)$ works very well, without any need for manual adjustment. Overall, the proposed temporal artificial stress diffusion can be recommended as a simple upgrade and extension of existing solvers where it can improve their stability and robustness.