1 Introduction

The Oldroyd-B model is one of the classical rheological models describing certain class of viscoelastic fluids found in many practical applications including engineering, environmental or biomedical sciences. Although other models are often used, this one remains a kind of a prototype, sharing the same structure and some behavior with many more advanced and sophisticated models. One of the well known disadvantages of Oldroyd-B type (and similar) models is that numerical simulations of the fluids flows governed by these models become unstable in certain flow regimes, namely in the case of high Weissenberg numbers [3, 30]. This problem is known as the High Weissenberg Number Problem (HWNP) [23].

This numerical issue is treated in a different way by different authors, leading to number of specific stabilization and discretization techniques trying to overcome the limitations arising from HWNP, see for example [12, 19, 23, 32, 42, 43] or the review paper [1] containing many additional references. One of the most efficient methods to overcome the HWNP is based on the logarithmic transformation for the conformation tensor [13], which allows to numerically capture the exponential growth of tensor components in certain regions. This method was widely used and refined by many authors [10, 11, 14, 24, 40]. This approach however requires a specialized code to be developed to implement the complete algorithm.

This is why certain simple stabilization algorithms are still of interest, due to their easy implementation into existing general CFD codes. One of the simplest and most widely used stabilization methods (not just for viscoelastic fluids flows) is based on the introduction of the so-called numerical viscosity (or diffusion) [16, 21, 22, 45]. In this approach an extra (artificial) numerical diffusion is added to the model at the discretization stage, leading to more stable simulations, i.e. more robust numerical solvers. However, this term is of purely numerical origin and has no physical significance [7].

The extra numerical diffusion can either be added as a separate term (to be discretized) within the model equations [20, 25] or it is naturally present while using low order discretization methods (e.g. as upwinding) [18, 26, 27]. This approach was successfully used during past decades in many numerical fluids flows simulations, for example the solution of hyperbolic problems including the Euler equations in gas dynamics [21, 28, 29, 44].

Due to simplicity and widespread use of the numerical diffusion technique, it is natural to try to use it also for treating the numerical instability of the Oldroyd-B fluids flows. In viscoelastic model however, the main problem is in the coupled system of transport-like equations for components of the viscoelastic stress tensor. The added numerical diffusion should thus be a tensorial quantity. The added diffusive term is typically directly proportional to the Laplacian of the extra stress tensor, i.e. \(\alpha \cdot \Delta \varvec{\tau }\) (using the notation introduced below). This approach is rather elementary, but it has some advantages. Most importantly, there exist some physical (or thermodynamical) arguments that a diffusive term of this type should be present in the model [2], which leads to a so called diffusive Oldroyd-B model [8]. So, in such case the added stress diffusion is physical feature and the diffusion coefficient \(\alpha\) is physical quantity coming from experimental measurements. So in fact this added tensorial diffusion leads to another rheological model and as a side effect, the increased stability of the model numerical solution is gained. This also indicates the drawback of using the added stress diffusion as a stabilization technique for the standard (non-diffusive) Oldroyd-B fluid flow simulations. The solution obtained using the extra diffusion term leads to different solution, corresponding to the diffusive Oldroyd-B model, rather than to the original (non-diffusive) Oldroyd-B model for which a solution was sought. In other words, the increased stability of the numerical solver (using \(\alpha > 0\)) comes at the expense of lower accuracy of the solution of the original Oldroyd-B model problem. The higher values of \(\alpha\), the further from the solution of the non-diffusive Oldroyd-B model are the numerical solutions obtained from of (other) diffusive Oldroyd-B model.

This means that in case the added diffusion term is used just to stabilize the standard (non-diffusive) Oldroyd-B model, the added diffusion is purely artificial and thus the coefficient \(\alpha >0\) should be kept as small as possible for the solution to remain as close as possible to the solution of the original non-diffusive problem. Technically speaking, when trying to solve Oldroyd-B flows at higher Weissenberg numbers, the numerical diffusion coefficient \(\alpha\) has to be increased, which will (soon or later) lead to a solution which is unacceptably far from the solution of the non-diffusive Oldroyd-B model.

Our previous experience with the simulation of Oldroyd-B like fluids flows indicate that the numerical instability develops during the initial stage of iterative process (towards steady state), which leads to the need of stronger stabilization during this initial phase of numerical simulation. It was also found that any level of added stress diffusion visibly affects not just the numerical stability of the solver, but also the final numerical solution itself. The solution of the model with added diffusion is typically smoother, with lower gradients and cut-off peaks of flow quantities (not just stress tensor) [8] or Sect. 5.1 in this paper. These initial observations led us to the idea of using modified stress diffusion term, that will mainly be active during initial steps of the iteration process, while it will be decaying later, potentially completely diminishing in the steady state.

We have proposed and tested several variants of the vanishing artificial stress diffusion. First variant was an added term proportional to Laplacian of time derivative of the stress tensor, i.e. term of the form \(\alpha \cdot \Delta \varvec{\tau }_{t}\). The logic behind this choice follows the original goal to have a stabilization term that will mainly be active at the beginning of the iterative process, when the pseudo-time derivatives of the solution are higher, while it automatically vanishes when the steady solution is reached. This concept proved to be quite successful, as we have shown in our previous work [34]. As one of the drawbacks of this kind of added stabilization term might be seen the fact that this extra term has a non-standard and non-physical form containing third order mixed partial derivative of the computed (tensorial) quantity. This might complicate some theoretical considerations justification of this model as well as its numerical implementation.

Therefore we were looking for some alternative stabilization terms that will have similar behavior, but can better be supported by mathematical theory [5, 9] and will be easier to implement into existing numerical solvers. Two different alternatives are presented and discussed in this paper. Both of them use the added Laplacian of the extra stress tensor multiplied by a diffusion coefficient \(\alpha\), which is no more a constant, but is constructed as a suitable function of time or time-derivative of the solution, i.e. the added terms have the form \(\alpha (t)\!\cdot \!\Delta \varvec{\tau }\) or \(\alpha (\partial _{t})\!\cdot \!\Delta \varvec{\tau }\). The results obtained using these two basic variants of the stabilization term are presented and discussed in detail. The focus is on their comparative assessment based on the qualitative behavior of the solution obtained by a Finite Element Method (FEM).

This method is one of the most commonly used in numerically solving problems described by partial differential equations [37, 46]. The numerical solution is constructed from a system of equations reformulated in the variational form. At discrete level it leads to decomposition of the physical domain of interest into a finite number of simple subdomains - the finite elements. On each element a set of nodes and low degree polynomials is defined, allowing to construct finite-dimensional function spaces for numerical approximation of the discretized problem. During the past decades the FEM became also a standard tool in investigation and numerical modeling of viscoelastis fluids flows problems [15, 31, 38, 39, 41].

This paper is continuation and extension of our previous works [4, 33, 34], therefore in order to minimize possible repetition of previously published information, we skip the detailed description of the weak formulation of the problem, numerical discretization and solver details which can be found in and [33,34,35]. Here we only review some essential details of the model and focus on the description of certain types of modified artificial stress diffusion terms.

The aim of this paper is to find out whether for given rheological model, numerical method and computational grid, the added artificial tensorial diffusion can help to shift (increase) the critical Weissenberg number. This is why in this paper we don’t explore other effects related to the use of alternative rheological models, numerical methods or grids. Here the focus is purely on the effects of vanishing the added artificial tensorial diffusion.

The structure of the paper will be as follows: After the Introduction, the mathematical model for Oldroyd-B fluids flow simulations is summarized in the Sect. 2, where also the references for numerical solver and its implementation are given. The different types of artificial stress diffusion are presented in Sect. 3, pointing out some main features and expected behavior. The whole Sect. 4 is dedicated to description of the solved test case and mainly to presentation of the obtained results and their assessment. The Sect. 5 concludes the paper by summarizing our experience coming from numerical results and outlining some possible future work in this area.

2 Mathematical model and numerical method

Mathematical description of flow of incompressible Oldroyd-B type fluid in two- or three- dimensional space can be written in the following non-dimensional form [3, 38, 39] or [41].

$$\begin{aligned} \left\{ \begin{array}{l} \mathcal {R}e \left( \dfrac{\partial \varvec{u}}{\partial t}+ \varvec{u}\cdot \nabla \varvec{u}\right) + \nabla p = 2(1-\eta )\nabla \cdot \varvec{\textsf{D}}+ \nabla \cdot \varvec{\tau }+\varvec{f}\\ \nabla \cdot \varvec{u}=0\\ \varvec{\tau }+\mathcal {W}e\left( \dfrac{\partial {\varvec{\tau }} }{\partial t}+{ \varvec{u}}\cdot \nabla \varvec{\tau } -\nabla \varvec{u}^T \cdot \varvec{\tau }- \varvec{\tau }\cdot \nabla {\varvec{u}}\right) =2\eta \varvec{\textsf{D}}\end{array} \right. \end{aligned}$$
(1)

where \(\mathcal {R}e = \dfrac{UL}{\nu }\) and \(\mathcal {W}e = \dfrac{\lambda U}{L}\) are the Reynolds and Weissenberg numbers, respectively. The viscoelastic extra stress (symmetric) tensor is denoted by \(\varvec{\tau }\). It contributes to the stress tensor \(\varvec{\textsf{T}}=2(1-\eta )\varvec{\textsf{D}}+ \varvec{\tau }\), where \(\varvec{\textsf{D}}=\dfrac{1}{2}\left( \nabla u + \nabla u^T\right)\) is the rate of strain tensor (symmetric part of velocity gradient). The external body force is denoted by \(\varvec{f}\).

The governing system is considered (and solved) in its unsteady form as the time-marching procedure is employed to search for steady solution as a limit of the unsteady iterations for pseudo-time \(t\longrightarrow \infty\). This approach allows us to use the vanishing-in-time stabilization terms.

The model is solved numerically using the so-called reference viscosity scheme. The extra stress tensor is split according to \(\varvec{\textsf{T}}=2(1-\eta )\varvec{\textsf{D}}+ \varvec{\tau }\) into the viscoelastic stress \(\varvec{\tau }\) and the viscous part \(\varvec{\tau }_s= 2(1-\eta )\varvec{\textsf{D}}\). The parameter \(\eta \in [0,1]\) is, in this case, the dimensionless polymer viscosity contribution. This splitting allows to decouple the kinematic and the non-Newtonian viscoelastic stress where the divergence of \(\varvec{\tau }\) is included in the momentum equations as a pseudo-body force while the constitutive equation includes a contribution from the Newtonian part.

More details on numerical methods and their implementation can be found in our previous works [33, 34], Chapter [35] in [6] or in the documentation of the solver package used [17].

3 Stabilization terms

As mentioned in the Introduction, several diffusion stabilization strategies were proposed and tested [8, 15]. One of them is based on the introduction of the so-called numerical viscosity (or diffusivity) [45] which is of purely numerical origin and has no physical meaning [7]. Following the ideas from [16, 21, 22], an extra (artificial) numerical diffusion term in the form \(\alpha \cdot \Delta \varvec{\tau }\) is added at the discretization stage to the transport-type equation for viscoelastic stress tensor [20, 25], resulting to a more stable numerical model. This additional tensorial diffusion leads to a model similar to the diffusive Oldroyd-B model [8] whose rheological behavior is different. It means that the solution obtained by the stabilized model using the extra diffusion term provides in fact a solution corresponding rather to the diffusive Oldroyd-B model, instead of the original (non-diffusive) model, for which the solution was sought. This is why in the stabilization of the standard (non-diffusive) Oldroyd-B model, the added diffusive term is purely artificial and therefore the coefficient \(\alpha > 0\) should always be kept as small as possible, in order to guarantee that the obtained (stabilized) solution remains close to the solution of the original non-diffusive problem. It should be noted that the need of stabilization only emerges when solving the Oldroyd-B fluid flows at larger Weissenberg numbers, meaning that the numerical diffusion coefficient \(\alpha\) must be increased in that regime. However this increase will for high Weissenberg numbers lead (sooner or later) to a solution that is unacceptably far from the solution of the non-diffusive Oldroyd-B model.

In one of our previous works [33] we have proposed and successfully tested an added numerical diffusion term proportional to Laplacian of time-derivative of the stress tensor in the form \(\alpha \cdot \Delta \varvec{\tau }_t\), where the coefficient \(\alpha\) was a suitably chosen constant. At discrete level the term was implemented as \(\alpha \cdot \Delta \left( {\varvec{\tau }}^n_h - {\varvec{\tau }}^{n-1}_h\right)\), meaning that Laplacian of the difference between two successive iterations of the discrete solution (at time levels n and \((n-1)\)) was used. It is easy to see that this added term vanishes in the limit steady case when \(\left( {\varvec{\tau }}^n_h - {\varvec{\tau }}^{n-1}_h\right) \rightarrow 0\), i.e. \(\varvec{\tau }^n_ h\rightarrow \varvec{\tau }\). So the added term is only present during the time marching procedure, but disappears once the solution converges to the steady state.

The numerical tests showed that this added diffusive term stabilizes the simulations for wide range of choices of the parameter \(\alpha\). Evidently the obtained steady solution corresponds to the solution of the original problem without diffusion. This is in contrast with the use of the standard diffusion term \(\alpha \cdot \Delta \varvec{\tau }\), where the added diffusive term doesn’t depends on time (or time-derivative) and therefore doesn’t in general vanishes at the steady state and thus the obtained solution rather corresponds to another, diffusive model [8]. This situation was demonstrated and discussed in our recently published paper [36]. Such additional stabilization term can possibly be interpreted as a kind of a residual smoothing technique, where the expression \(\left( {\varvec{\tau }}^n_h - {\varvec{\tau }}^{n-1}_h\right)\) represents a steady residual of the problem solved.

The above mentioned stabilization term worked well, but keeping in mind some of the open theoretical questions related to well posedness of the modified model, we kept looking for some alternative stabilization terms that will be closer to the standard numerical diffusion \(\alpha \cdot \Delta \varvec{\tau }\). For such simple diffusive model some theoretical results are available discussing the existence, uniqueness or regularity of the solution [9].

This is why we have focused our investigation on simple generalization of this standard diffusive model, by making the diffusion coefficient \(\alpha (\cdot )\) variable, depending either just on time or on some time-derivative of the solution.

The constitutive equation including such generalized diffusion can be written in the form

$$\begin{aligned} \varvec{\tau }+\mathcal {W}e\left( \dfrac{\partial {\varvec{\tau }} }{\partial t}+{ \varvec{u}}\cdot \nabla \varvec{\tau } -\nabla \varvec{u}^T \cdot \varvec{\tau }- \varvec{\tau }\cdot \nabla {\varvec{u}}\right) =2\eta \varvec{\textsf{D}}+\alpha (\cdot )\Delta \varvec{\tau } \end{aligned}$$
(2)

for which the weak formulation given by

$$\begin{aligned}{} & {} \displaystyle \int _\Omega \left[ \varvec{\tau }+\mathcal {W}e\left( \dfrac{\partial {\varvec{\tau }} }{\partial t}+{ \varvec{u}}\cdot \nabla \varvec{\tau } \right) \right]:\textbf{S}= \\{} & {} \quad = \displaystyle \int _\Omega \! \left[ 2\eta \varvec{\textsf{D}}+\mathcal {W}e\left( \nabla \varvec{u}^T \cdot \varvec{\tau }+ \varvec{\tau }\cdot \nabla {\varvec{u}}\right) \right] :\textbf{S} - \alpha (\cdot ) \nabla \varvec{\tau }:\nabla \mathcal {S},{} & {} \quad\forall \textbf{S}\! \in \! \mathcal {S} \end{aligned}$$

where \({\mathcal {S}} =\left\{ \textbf{S}\in [L^2(\Omega )]^{2\times 2}: \textbf{S}^T=\textbf{S}\right\}\).

The approximate finite element problem for the tensor can be reformulated as:

For each \(t\in [0,T_{ f }], \ \varvec{\tau }_{0,h} \in S_h = \left[ \mathcal {M}_h\right] ^{2\times 2}\), find \(\varvec{\tau }_{h} \in \mathcal {S}_h\) such that each component \(\varvec{\tau }_{h,ij}\) verifies

$$\begin{aligned}{} & {} \displaystyle \int _\Omega \left(\varvec{\tau }^n_{h,ij} + \mathcal {W}e\dfrac{{\varvec{\tau }}^n_{h,ij} - {\varvec{\tau }}^{n-1}_{\star ,h,ij} }{\Delta t}\right) :\textbf{S}_{h,ij}= \nonumber \\{} & {} \quad \displaystyle =\int _\Omega \!\! \left[ 2\eta \varvec{\textsf{D}}\!+\!\mathcal {W}e\left( \nabla \varvec{u}^T \!\!\cdot \varvec{\tau }^n_h\!\!+\!\! \varvec{\tau }^n_h\!\cdot \!\nabla {\varvec{u}}\right) \right] _{ij}\!\!: \textbf{S}_{h,ij} - \alpha (\cdot ) \nabla \varvec{\tau } ^n_{h,ij} \! \!:\!\textbf{S}_{h,ij} \end{aligned}$$
(3)

for all \(\textbf{S}_h\in \mathcal {S}_h\), where \({\tau }^{n-1}_{\star ,h} = {\tau }^{n-1}_h(x^\star )\), being \(x^\star\) the position at time \(t_{n-1}\) of the fluid parcel situated at x at time \(t_n\) (Characteristic Galerkin Method [33, 46]). The same kind of discretization was also used for momentum equations. The complete variational formulation of the model and its finite element discretization was described in [35].

The vanishing diffusion model doesn’t changes significantly with respect to the original non-diffusive model, and also the numerical implementation remains almost identical, except the coefficient \(\alpha (\cdot )\) now being a suitably chosen function. In this paper two different types of vanishing artificial diffusion coefficients are described and tested:

  1. 1.

    Time dependent function \(\alpha (t)\) – Diffusion coefficient \(\alpha\) is monotonically decaying with (pseudo) time t and vanishing for the limit case of \(t\longrightarrow \infty\), where the steady state solution should be reached. Different shapes of such monotonically decaying functions were considered and tested, searching for optimal initial values of the diffusion coefficient and suitable decay rate. The function \(\alpha (t)\) used in the presented simulation has the form

    $$\begin{aligned} \alpha (t)=\alpha _{0}\cdot \frac{1}{1+\varepsilon \cdot t}\quad , \end{aligned}$$
    (4)

    where \(\alpha _{0}=\alpha (0)\) is the initial value for the diffusion parameter \(\alpha\), while the adjustable parameter \(\varepsilon\) affects the rate of decay of the function (by scaling the time variable).

  2. 2.

    Time-derivative dependent function \(\alpha (\phi _t)\) – Diffusion coefficient \(\alpha\) is made proportional to (dependent on) the time-derivative of some solved quantity \(\phi\). The dependence is such that the diffusion coefficient \(\alpha \longrightarrow 0\) for the steady state solution where \(\phi _{t}\longrightarrow 0\). Proposed and tested were different functional dependencies of \(\alpha\) on the time derivative \(\phi _{t}\) as well as different choices of the indicator variable \(\phi\) for evaluation of the time derivative (e.g. pressure p, tensor components \(\tau _{ij}\) or tensor norm \(\Vert \varvec{\tau } \Vert\)). The general form of the functional dependency of \(\alpha (\phi _t)\) is the following:

    $$\begin{aligned} \alpha (\phi _t)=\alpha _{0}\cdot \frac{\varepsilon }{\varepsilon +(1-\varepsilon )\cdot \left( \Vert \phi _{t}\Vert \right) ^{m}}, \end{aligned}$$
    (5)

    where again \(\alpha _{0}=\alpha (0)\) is the initial value for the diffusion parameter \(\alpha\), and the adjustable parameters \(\varepsilon\) and m affect the rate of decay of the function.

The exact form of the functional dependencies (4) and (5) was chosen ad-hoc, only keeping in mind suitable (adjustable) decay rate and limit behavior for \(t\rightarrow \infty\) and \(\phi _{t}\rightarrow 0\) respectively. Other alternative functions we have used led to similar results and had no significant impact on the presented findings and conclusions.

4 Numerical setup

The described mathematical model including the added stabilization terms was solved using the FreeFem++ [17] toolboxes. For all simulations the Reynolds number was fixed to \(\mathcal {R}e= 1000\), the dimensionless elastic viscosity parameter \(\eta = 0.1\). The Weissenberg number \(\mathcal {W}e\) varies from zero (for Newtonian fluid) up to certain critical value, for which the solution becomes unstable and solver fails. To obtain results for higher Weissenberg numbers, the continuation method is applied, meaning that the Weissenberg number is increased incrementally, by taking the converged solution obtained for lower \(\mathcal {W}e\) as an initial condition for the simulation at higher \(\mathcal {W}e\).

4.1 Computational domain

Numerical simulations focus on a 2D corrugated channel test case we have previously used for evaluation of other stabilization strategies [4, 33, 34]. The corrugated channel (pipe) consists of (sufficiently long) straight inlet and outlet parts with dimension (width/height) D. In the middle of the channel length are smoothly attached several identical sinusoidally shaped segments. Dimensions of these segments are defined by their minimum and maximum diameters \(D_{min}\) resp. \(D_{max}\) and the segment length \(L_{seg}\) (Fig. 1).

Fig. 1
figure 1

Geometry definition for the corrugated pipe case

Full size image

For this domain a quasi regular unstructured grid (triangulation) was generated by the Delaunay-Voronoi algorithm [7], using at the boundary 10 points per unit length (i.e. per D) resulting into a mesh of diameter around 0.0844 consisting of 3070 triangular elements, leading to 6481 nodes for \(P_2\) elements and 1706 nodes for \(P_1\) elements.

The results presented in this paper were obtained on single, identical grid used in all simulations. This was done intentionally, trying to find out if for given computational case and code, the range of attainable Weissenberg numbers can be extended by simply adding an artificial diffusion term.

4.2 Boundary conditions

The boundary conditions for the considered case are quite standard. The fully developed profile for velocity and extra stress tensor is used at the inlet, imposed as Dirichlet condition using analytical solution of Poisseuille flow of the Oldroyd-B fluid obtained for given \(\mathcal {R}e\) and \(\mathcal {W}e\). No-slip conditions are used for velocity on the channel walls. No boundary conditions are explicitly imposed on the extra stress tensor components. In case of the added stress tensor (artificial) diffusion, the finite element solution naturally satisfies the homogeneous Neumann conditions for the stress tensor. The steady solution was sought using the time-marching procedure, i.e. by considering the steady (fixed in time) boundary conditions, and searching for the limit steady solution for \(t\rightarrow \infty\). No diffusion specific modifications of the boundary conditions were made.

4.3 Artificial diffusion setup

For both the time dependent and time-derivative dependent functions \(\alpha (\cdot )\) the same initial value \(\alpha (0) =\alpha _{0}= 10^{-4}\). This allows for mutual comparison of the two stabilization approaches as well as the comparison with a reference solution obtained using the constant diffusion coefficient \(\alpha =\alpha _{0}\).

For the time dependent artificial diffusion coefficient \(\alpha (t)\) in the form (4), different values of the decay parameter \(\varepsilon\) was tested in the range from \(\varepsilon =0\) (for constant \(\alpha\)) up to \(\varepsilon =1.0\) when the function \(\alpha (t)\) decays fastest (Fig. 2).

Fig. 2
figure 2

Behavior of the function \(\alpha (t)\) (diffusion coefficient) with respect to time

Full size image

In case of the time-derivative dependent diffusion coefficient \(\alpha (\phi _{t})\) the decay parameter \(\varepsilon\) was set to 0.01, while the power \(m=5\). As indicator variable \(\phi\) (from which the time derivative is evaluated), we have chosen and tested the pressure p, extra stress tensor components \(\tau _{ij}\) and the norm of the extra stress tensor \(\Vert \varvec{\tau } \Vert\).

5 Numerical results

The flow in the above described two-dimensional channel was solved to assess some of the essential performance and robustness characteristics associated with the proposed stabilization methods. The results of simulations are split into two separate sections, first for the time dependent and second for the time-derivative dependent stabilization diffusion coefficients. Separate discussion of results is presented within these two sections, while the general assessment comments are left for the Conclusions section at the end of the paper.

The results are presented in the form of graphs of selected quantities (stress components, stress tension) along the wall. These quantities are of practical importance, but also the solution gradients are typically highest close to the wall and thus the solution instabilities as well as the solution smoothing is most apparent there. In addition to these graphs we provide some plots of the norms of the Laplacian of the stress tensor (documenting the solution smoothness) and the levels of maximum stress for different cases and stabilization parameters. To assess the influence of the added (non-physical/artificial) diffusion term \(\alpha \Delta {\varvec{\tau }}\), its norm is also evaluated and presented for individual cases.

5.1 Time dependent diffusion coefficient

Many results for this case were already presented in our previous work [34], so only a selection of few results is shown here to support the discussion presented in this paper.

Figure 3 shows how the evolution (in iterative time) of the extra stress tensor field is affected by the different choices in the decay rate of the diffusion coefficient \(\alpha (t)\). The norm of the Laplacian \(\Vert \Delta \varvec{\tau }\Vert\) of the extra stress tensor is considered as a measure of the smoothness of the tensor field. The comparison of results for low Weissenberg number \(\mathcal {W}e = 0.1\) (on the left) and close to critical \(\mathcal {W}e = 0.4\) (on the right) shows significant increase (by an order of magnitude) of the values of \(\Vert \Delta \varvec{\tau }\Vert\), which demonstrates how much different (and difficult) are the simulations for low versus high Weissenberg numbers. The high Weissenberg number case is substantially “less smooth” and thus more difficult to approximate numerically.

For the low Weissenberg number \(\mathcal {W}e = 0.1\) the solution can be obtained without any added stabilization, which corresponds to setting \(\alpha =0\). Such solution (of the original, non-diffusive problem) can be used for comparison with solutions obtained using different artificial diffusion settings. The values obtained from the non-diffusive case are those we want to obtain with artificial stabilization.

Fig. 3
figure 3

Comparison of the \(\Vert \Delta \varvec{\tau }^n\Vert\) for different values of parameter \(\varepsilon\) of diffusion coefficient \(\alpha (t)=\frac{10^{-4}}{1+\varepsilon t}\) for the Weissenberg number \(\mathcal {W}e = 0.1\) (left) and \(\mathcal {W}e = 0.4\)(right)

Full size image

The vanishing diffusion concept is based on idea that the non-physical added numerical diffusion term \(\alpha (\cdot )\Delta \varvec{\tau }\) vanishes for \(t\rightarrow \infty\). This assumption is verified in the Fig. 4 showing the norm of the whole artificial diffusion term \(\Vert \alpha (t)\Delta \varvec{\tau }^n\Vert\) during the iterative process. It clearly shows that the diffusion coefficient function \(\alpha (t)\) decay rate parameter \(\varepsilon\) should be set high enough to guarantee that not only \(\alpha (t)\) will decay, but also the whole artificial diffusion term will decay in time and vanish for \(t\rightarrow \infty\) (will be negligible at the moment the iterations are stopped).

Fig. 4
figure 4

Comparison of the \(\Vert \alpha (t)\Delta \varvec{\tau }^n\Vert\) for different values of \(\varepsilon\) in the diffusion coefficient \(\alpha (t) = \frac{10^{-4}}{1+\varepsilon t}\) for the Weissenberg number \(\mathcal {W}e = 0.1\) (left) and \(\mathcal {W}e = 0.4\)(right)

Full size image

Figure 5 is showing the final values of the \(\Vert \Delta \tau ^n\Vert\) and \(\Vert \alpha (t)\Delta \tau ^n\Vert\) for a range of Weissenberg numbers \(\mathcal {W}e\). Again the \(\Vert \Delta \tau ^n\Vert\) represents the solution (non-)smoothness and \(\Vert \alpha (t)\Delta \tau ^n\Vert\) shows whether the non-physical diffusion term vanished. The horizontal extent of individual graphs in Fig. 5 shows what is the maximum attainable (critical) Weissenberg number for the given artificial diffusion setting. Evidently the solutions with non-vanished artificial diffusion exist for higher Weissenberg numbers (than without the stress diffusion), however they are much smoother (with lower \(\Vert \Delta \tau ^n\Vert\)) than the reference non-diffusive solution of the original problem.

Fig. 5
figure 5

Comparison of the \(\Vert \Delta \varvec{\tau }^n\Vert\) (left) and \(\Vert \alpha (\phi _t)\Delta \varvec{\tau }^n\Vert\) (right) for different values of parameter \(\varepsilon\) of diffusion coefficient \(\varepsilon\) of \(\alpha (t)=\frac{10^{-4}}{1+\varepsilon t}\) depending on Weissenberg numbers \(\mathcal {W}e\)

Full size image

More detailed information about the local smoothing properties of the tensorial stress diffusion can be found in the Fig. 6 showing the plots of individual stress tensor components along the channel wall. It is clearly visible that while the smooth part of the graphs is well resolved for all settings the curve peaks (local minima and maxima) are heavily affected by the applied numerical diffusion. The high constant or not yet vanished time-dependent diffusion coefficients \(\alpha\) lead to oversmoothed solution, where the extremal values clearly differ from those found in the reference non-diffusive solution (obtained for \(\alpha (t)=0\)).

Fig. 6
figure 6

Behavior of the extra stress component \(\tau\) along the channel wall for different values of parameter \(\varepsilon\) of diffusion coefficient \(\alpha (t)=\frac{10^{-4}}{1+\varepsilon t}\) for the Weissenberg numbers \(\mathcal {W}e = 0.1\) (left) and \(\mathcal {W}e = 0.4\) (right)

Full size image

The values of maxima and minima of tensor components on the channel wall are shown in the Figs. 7 and 8 respectively. The maximum and minimum values of the original problem solution (without artificial diffusion) are only preserved for the fastest decaying diffusion coefficient \(\alpha (t)\). For all the other settings the solution is significantly altered.

Fig. 7
figure 7

Maximum values of the extra stress tensor components along the channel wall for different setting of parameter \(\varepsilon\) in the diffusion coefficient \(\alpha (t)=\frac{10^{-4}}{1+\varepsilon t}\), depending on the Weissenberg number \(\mathcal {W}e\)

Full size image
Fig. 8
figure 8

Minimum values of the extra stress tensor along the channel wall for different setting of the parameter \(\varepsilon\) in the diffusion coefficient \(\alpha (t)=\frac{10^{-4}}{1+\varepsilon t}\), depending on the Weissenberg number \(\mathcal {W}e\)

Full size image

From the physical (technical) point of view an important quantity to look for is the stress tension on the wall \(-(\boldsymbol{\tau} \cdot \bf{n})\cdot \bf{t}|_w\). This quantity is plotted in the Fig. 9 with respect to axial coordinate, for different settings of the decay parameter \(\varepsilon\) in the diffusion coefficient \(\alpha (t)\). Again it is clear that the peaks in the stress tension are unresolved for slowly decaying (and thus not yet vanished) stress diffusion coefficients \(\alpha (t)\). Such underestimation of the forces acting on the channel wall may have significant impact in technical applications.

Fig. 9
figure 9

Stress tension on the wall \(-(\boldsymbol{\tau} \cdot \bf{n})\cdot \bf{t} |_w\) along the channel wall for different values of parameter \(\varepsilon\) in the diffusion coefficient \(\alpha (t)=\frac{10^{-4}}{1+\varepsilon t}\), for the Weissenberg numbers \(\mathcal {W}e = 0.1\) (left) and \(\mathcal {W}e = 0.4\) (right)

Full size image
Fig. 10
figure 10

Plots of the maximum (left) and minimum (right) values of stress tension \(-(\boldsymbol{\tau} \cdot \bf{n})\cdot \bf{t}|_w\) along the channel wall for different values of parameter \(\varepsilon\) in the diffusion coefficient \(\alpha (t)=\frac{10^{-4}}{1+\varepsilon t}\), depending on the Weissenberg number \(\mathcal {W}e\)

Full size image

The underestimation of the extremal values of the stress tension on the wall was confirmed for the whole range of Weissenberg numbers, as it can be seen from the Fig. 10

Another quantity of technical interest for the considered cases is the pressure drop between the inlet and outlet of the channel. It corresponds to force (or power) needed to push the fluid through the channel. This information is crucial in technical as well as in biomedical applications of the model. The pressure drop versus Weissenberg number plots are shown in the Fig. 11 for various setting of the diffusion coefficients \(\alpha (t)\). Evidently the artificial diffusion affects the predicted pressure drop values and only the fastest decay of \(\alpha (t)\) provides the desired solution comparable with the original problem without the artificial diffusion.

Fig. 11
figure 11

Comparison of the pressure drop for different values of parameter \(\varepsilon\) of diffusion coefficient \(\alpha (t)=\frac{10^{-4}}{1+\varepsilon t}\) depending on Weissenberg numbers \(\mathcal {W}e\)

Full size image

The quantitative differences between the non-diffusive model solution and the solution using the artificial diffusive term can be assessed by computing the norm of differences between the two solutions (which can be considered as an error with respect to Oldroyd-B model). From the quantities obtained in simulations we present in the Fig. 12 the differences of the shear component \(\tau _{12}\) of the stress tensor.

Fig. 12
figure 12

Norm of the difference of shear stress component \(\tau _{12}\) between the diffusive and non-diffusive model solution, depending on the Weissenberg number \(\mathcal {W}e\), for different values of the decay parameter \(\varepsilon\) in the artificial diffusion coefficient \(\alpha (t)=\frac{10^{-4}}{1+\varepsilon t}\)

Full size image

Evidently just for the fastest decay (\(\varepsilon =1.0\)) the results of the original non-diffusive model are comparable with the results obtained using the artificial diffusion term. This is confirmed in the Table 1, where the numerical values corresponding to Fig. 12 are shown.

Table 1 The shear stress differences \(\left\| \tau _{12,D}- \tau _{12,N}\right\| _2\) between the diffusive and non-diffusive model, depending on \(\mathcal {W}e\) for various parameters \(\varepsilon\) in the artificial diffusion coefficient \(\alpha (t)\)
Full size table

The same conclusion can also be drawn when comparing the pressure drop for the original non-diffusive Oldroyd-B model with the results obtained using the artificial diffusion term. Figure 13 shows that the difference is bigger when the decay rate parameter \(\varepsilon\) is decreased. Moreover such solution differences are getting more pronounced for higher Weissenberg numbers, where the additional artificial diffusion is supposed to be used.

Fig. 13
figure 13

Pressure drop difference between the diffusive and non-diffusive model solution, depending on the Weissenberg number \(\mathcal {W}e\), for different values of the decay parameter \(\varepsilon\) in the artificial diffusion coefficient \(\alpha (t)=\frac{10^{-4}}{1+\varepsilon t}\)

Full size image

The numerical values of the pressure drop differences (errors) shown in the Fig. 13 are summarized in the Table 2.

Table 2 The value of pressure drop differences \(\left| \Delta P_D- \Delta P_N\right|\) between the non-diffusive and diffusive model, depending on \(\mathcal {W}e\) for various parameters \(\varepsilon\) in the artificial diffusion coefficient \(\alpha (t)\)
Full size table

5.2 Time-derivative dependent diffusion coefficient

The time-derivative dependent artificial diffusion coefficient is for the first time tested here, so the core of the results presented in this paper is dedicated to this specific case. The preliminary simulations we have performed using the artificial diffusion with the variable coefficient \(\alpha (\phi _t)\) of the form (5) have shown that the performance of this kind of stabilization is strongly affected by two factors.

First, it is crucial to decide which quantity \(\phi\) is used as an indicator to evaluate the time derivative for calculating the diffusion coefficient \(\alpha\). From the variants we have tested, the results for \(\phi =p\) (pressure), and the norm of the extra stress tensor \(\phi =\Vert \varvec{\tau } \Vert\) are presented and compared here.

Second, it became apparent that the use of instantaneous values of time-derivatives just from the previous iteration can lead to destabilization of the code. A much better choice seems to be to consider some recent history of \(\alpha (\phi _t)\), in the form of a floating average, to establish the actual value of \(\alpha\) for the next iteration. Therefore the tests were performed for the \(\alpha\) being a result of averaging over the previous L iterations, starting from \(L=10\) and increasing the length of the averaging history always by factor 2, leading to series of tests for \(L=10,\,20,\,40,\,80,\,160,\,320,\,640,\,1280\). The resulting averaged diffusion coefficient is thus defined by

$$\begin{aligned} \alpha (\phi _t)=\frac{1}{L}\sum _{k=n-L}^{n}\alpha (\phi _t^k)=\frac{\alpha _{0}}{L}\sum _{k=n-L}^{n}\frac{\varepsilon }{\varepsilon +(1-\varepsilon )\cdot \left( \Vert \phi _{t}^{k}\Vert \right) ^{m}}\quad , \end{aligned}$$
(6)

where k stands for the general time level index, n is the actual time level, L is the number of time levels over which the diffusion coefficient is averaged. As before the \(\alpha _0\), \(\varepsilon\) and m are parameters defining the shape of the function \(\alpha\).

The results presented hereafter are just selection from the series of tests performed, showing the major differences between the stabilization methods and effects of their parameters adjustments.

Figure 14 shows the pseudo-time history of comparison of the norm of the Laplacian of extra stress tensor \(\Vert \Delta \varvec{\tau } \Vert\) obtained using the stabilization based on the choice \(\phi =\Vert \varvec{\tau } \Vert\) (left column) and \(\phi =p\) (right column). Two observations can be made from this figure. First the solution obtained using \(\alpha (\Vert \varvec{\tau }_t\Vert )\) is less affected by the choice of the averaging history length L. Second, while the pressure based \(\alpha (\Vert p_t\Vert )\) is used, the Laplacian of the extra stress tensor is smaller (in the norm), which indicates smoother solution.

Fig. 14
figure 14

Comparison of \(\Vert \Delta \varvec{\tau }^n\Vert\) for diffusion coefficient \(\alpha (\Vert \varvec{\tau }_t\Vert )\) (left) and \(\alpha (\Vert p_t\Vert )\) (right) for different lengths of the averaging history L and range of Weissenberg numbers \(\mathcal {W}e\)

Full size image

The observations made in the Fig. 14 are just confirmed in the Fig. 15 showing the maximum of norm of the Laplacian of the extra stress tensor in dependence on the Weissenberg number \(\mathcal {W}e\). The stabilization using \(\alpha (\Vert \varvec{\tau }_t\Vert )\) (shown on left) performs well independently of the choice of the averaging history length L, all the curves are almost identical. On the other hand, while using the pressure based \(\alpha (\Vert p_t\Vert )\) (shown on right), the results heavily depend on the choice of L, leading to much lower maximum values of \(\Vert \Delta \varvec{\tau } \Vert\), which can indicate over-smoothed solution with extrema being significantly smeared. It is worth noting that some results were obtained for higher \(\mathcal {W}e\) using \(\alpha (\Vert p_t\Vert )\), while the stabilization using \(\alpha (\Vert \varvec{\tau }_t\Vert )\) already failed, but there is a doubt whether such over-smoothed results are still physically relevant to the solved case.

Fig. 15
figure 15

Comparison of the maximum of \(\Vert \Delta \varvec{\tau }\Vert\) for diffusion coefficient \(\alpha (\Vert \varvec{\tau }_t\Vert )\) (left) and \(\alpha (\Vert p_t\Vert )\) (right) depending on Weissenberg number \(\mathcal {W}e\) for different lengths of the averaging history L

Full size image

The trend is clear form Fig. 16 showing the pseudo-time dependence of the norm of the added artificial diffusion term, i.e. the \(\Vert \alpha (\phi _t)\Delta \varvec{\tau }^n\Vert\). It is evident that while the norm of the artificially added term really vanishes for the choice \(\alpha (\Vert \varvec{\tau }_t\Vert )\) (on the left), while for the pressure dependent \(\alpha (\Vert p_t\Vert )\) (on the right) the term is still visibly present and active even for very small Weissenberg numbers. This can probably by attributed to the fact that some temporal oscillations of pressure persist in the computational field, leading to non vanishing \(p_t\), which is even more pronounced when the averaging pseudo-time interval length L becomes larger.

Fig. 16
figure 16

Comparison of \(\Vert \alpha (\phi _t)\Delta \varvec{\tau }^n\Vert\) for diffusion coefficient \(\alpha (\Vert \varvec{\tau }_t\Vert )\) (left) and \(\alpha (\Vert p_t\Vert )\) (right) for different lengths of the averaging history L and range of Weissenberg numbers \(\mathcal {W}e\)

Full size image

Confirmation for this partial conclusion comes form the Fig. 17, showing the final value of the \(\Vert \alpha (\phi _t)\Delta \varvec{\tau }^n\Vert\) at the moment the simulation was stopped after finite number of iterations, for a range of Weissenberg numbers \(\mathcal {W}e\). Different lines in the figure correspond to different averaging history length L, while results for \(\alpha (\Vert \varvec{\tau }_t\Vert )\) and \(\alpha (\Vert p_t\Vert )\) are on the left and right side respectively.

Fig. 17
figure 17

Comparison of the \(\Vert \alpha (\phi _t)\Delta \varvec{\tau }^n\Vert\) for diffusion coefficient \(\alpha (\Vert \varvec{\tau }_t\Vert )\) (left) and \(\alpha (\Vert p_t\Vert )\) (right) depending on Weissenberg numbers \(\mathcal {W}e\) for different lengths of the averaging history L

Full size image

In Fig. 17, the behavior of the artificially added term is satisfactory for \(\alpha (\Vert \varvec{\tau }_t\Vert )\), the term remains small, even if some minor random looking oscillations appear in the graphs. However for the case with \(\alpha (\Vert p_t\Vert )\) the artificial term continues to be significantly large, which only gets worse for higher Weissenberg numbers and longer averaging history. To confirm whether the artificial diffusion coefficient \(\alpha (\phi _t)\) really vanishes in time, Fig. 18 shows the pseudo-time history of the \(\alpha (\Vert \varvec{\tau }_t\Vert )\) (on the left), and for the pressure dependent \(\alpha (\Vert p_t\Vert )\) (on the right). It clearly documents that the stress dependent diffusion coefficient small during the iterative process, vanishes for smaller Weissenberg numbers and its behavior is smoother for longer averaging history. On the other hand the pressure based \(\alpha (\Vert p_t\Vert )\) exhibits quite oscillatory behavior and moreover it remains significantly large when the history averaging is applied to it.

Fig. 18
figure 18

Pseudo-time evolution of \(\alpha (\phi _t)\) for diffusion coefficient \(\alpha (\Vert \varvec{\tau }_t\Vert )\) (left) and \(\alpha (\Vert p_t\Vert )\) (right) for different lengths of the averaging history L and range of Weissenberg numbers \(\mathcal {W}e\)

Full size image

So far we were just looking at some qualitative indicators of the solution, but the difference in effect of both versions of the time-derivative dependent stabilization term are also apparent on the final solution itself. Figure 19 shows the graphs of the numerically obtained stress component \(\tau _{11}\) along the channel wall. For convenience the channel wall geometry is shown at the bottom of each figure as well. The stress develops very sharp peaks along the wall, with higher values in higher Weissenberg number regimes. Again the results are compared for the two versions of the time-derivative dependent artificial stress diffusion (\(\alpha (\Vert \varvec{\tau }_t\Vert )\) on the left, \(\alpha (\Vert p_t\Vert )\) on the right), where both are used with the history averaging length \(L=20\), 80 and 320. For the stress dependent diffusion coefficient \(\alpha (\Vert \varvec{\tau }_t\Vert )\), the peaks are well resolved for all three values of L (even if some minor differences can be observed). But for the pressure based \(\alpha (\Vert p_t\Vert )\) leads to significant differences in the solution when the value of L is increased, leading to reduction of the extremal stress values on the wall and visibly (over)smoothed solution. This observation is in line with out previous conclusion concerning the non vanishing values of \(\alpha (\Vert p_t\Vert )\) and \(\Vert \alpha (\phi _t)\Delta \varvec{\tau }\Vert\) shown in Figs. 18 and 16 respectively. The same conclusion can be drawn from results for other lengths L and other two tensor components \(\tau _{12}\) and \(\tau _{22}\) (not shown here).

Fig. 19
figure 19

Behavior of the extra stress component \(\tau _{11}\) along the channel wall for diffusion coefficient \(\alpha (\Vert \varvec{\tau }_t\Vert )\) (left) and \(\alpha (\Vert p_t\Vert )\) (right) for different lengths of the averaging history L and range of Weissenberg numbers \(\mathcal {W}e\)

Full size image

In technical application typically rather than individual stress tensor components the stress tension (or wall shear stress) plays an important role. The corresponding graphs of the stress tension on the wall are shown in Fig. 20. Again the stronger smoothing of the stabilization term with \(\alpha (\Vert p_t\Vert )\) is apparent compared to the cases with \(\alpha (\Vert \varvec{\tau }_t\Vert )\). The performance of the pressure based stabilization becomes worse with longer history averaging L, but already for the shortest L the peaks in the solution are visibly smaller than in the results obtained with \(\alpha (\Vert \varvec{\tau }_t\Vert )\).

Fig. 20
figure 20

Behavior of the stress tension on the wall along the channel wall for diffusion coefficient \(\alpha (\Vert \varvec{\tau }_t\Vert )\) (left) and \(\alpha (\Vert p_t\Vert )\) (right) for different lengths of the averaging history L and range of Weissenberg numbers \(\mathcal {W}e\)

Full size image

This comparative advantage of the stress based artificial diffusion coefficient is confirmed when plotting the maximum values of the elastic stress tension along the wall for the two variants of the stabilization in Fig. 21.

Fig. 21
figure 21

Comparison of the maximum of the stress tension on the wall for diffusion coefficient \(\alpha (\Vert \varvec{\tau }_t\Vert )\) (left) and \(\alpha (\Vert p_t\Vert )\) (right) depending on Weissenberg numbers \(\mathcal {W}e\) for different lengths of the averaging history L.

Full size image

6 Conclusions and remarks

Numerical simulations were performed for two basic types of vanishing artificial stress diffusion under different parameter settings. In both types of diffusion, the standard Laplacian of stress tensor \(\Delta \varvec{\tau }\) was added to the model, multiplied by suitable diffusion coefficient \(\alpha\). In the first type the coefficient was considered as a function \(\alpha (t)\) explicitly depending on (pseudo-) time as in (4). The second type of vanishing artificial diffusion used the diffusion coefficient \(\alpha (\phi _t)\) defined as explicit function of the (pseudo-) time-derivative of certain indicator quantity as in (5). The conclusions arising from the obtained numerical results can thus be formulated separately for the two types of stabilization.

The time-dependent diffusion with \(\alpha (t)\) is the simplest method of implementation of vanishing artificial diffusion. The diffusion coefficient \(\alpha (t)\rightarrow 0\) as \(t\rightarrow \infty\) leading asymptotically to a “diffusion-free” model solution. In practical computations however the simulation is stopped at finite time t (after finite number n of iterations in (pseudo-) time), which means the decay rate of the function \(\alpha _t\) must be adjusted to guarantee small enough artificial diffusion at the moment the simulation is stopped and considered converged. This is not a problem for simulations with pre-specified fixed number of iterations, but it might be an issue, when the simulation stopping criterion is set for example based on the steady residual (proportional to pseudo time derivative of the solution) and we don’t know a-priori how fast the decay of \(\alpha (t)\) should be. In general the simulations with the stabilization using the variable diffusion coefficient \(\alpha (t)\) decaying monotonically in time have shown it’s efficiency when the initial value of \(\alpha\) and its decay rate were chosen properly, depending on the Weissenberg number and expected number of iterations to steady state. In such case the use of artificial tensorial diffusion was an efficient way of the stabilization, which was moreover safe to use, because this added artificial stabilization term vanished before the simulations stopped. This is definitely better behavior than the use of artificial diffusion stabilization with constant diffusion coefficient, where the final solution is always affected by the added term and different diffusion coefficients lead to different solutions. The main disadvantage of the first type of time dependent vanishing diffusion coefficient \(\alpha (t)\) is the pre-determined and fixed decay rate for \(\alpha\), just depending on time, without considering the actual solution or its steady residual. This is why we have proposed and intensively studied some other stabilizations depending on time derivative of the solution.

The time-derivative dependent diffusion with \(\alpha (\phi _{t})\) is the second, more complex type of stabilization we have used in the presented simulations. In this case the variable diffusion coefficient \(\alpha\) is made proportional to (norm of) time derivative of some indicator quantity \(\phi\) so that \(\alpha (\phi _t)\rightarrow 0\) for \(\phi _t\rightarrow 0\) (as \(t\rightarrow \infty\)). In this case the amount of added artificial diffusion is self-adjusted depending on the value of (pseudo-)time derivative of the solution and automatically vanishes when approaching the steady solution. This proved to be a nice and simple way of adjusting the artificial diffusion to the actual solution during the iterative process, without need to pre-define some solution independent decay rate for \(\alpha\). This automatic self-adjustment of of the diffusion coefficient \(\alpha\) helps namely in the case of higher Weissenberg numbers when the onset of temporal solution oscillations was often observed. This temporal oscillatory behavior leading to numerical instability can efficiently be suppressed by temporarily increasing the numerical diffusion. On the other hand, the automatic adjustment of the diffusion coefficient based on time-derivative rather than just on time brought some new problems. First, it wasn’t clear what quantity \(\phi\) should be used to indicate the need of adjustment of \(\alpha\) and be plugged in the formula for \(\alpha (\phi _t)\). After a number of preliminary simulations performed, we have decided to present results using two choices of the indicator variable \(\phi\). First is a flow field characteristic variable p corresponding to incompressible pressure and the second choice is the extra stress tensor norm \(\Vert \varvec{\tau } \Vert\). In both cases the chosen indicator variables are somehow related to forces (temporal force changes) acting on the fluid. The number of simulations presented in this paper for the case \(\phi =p\) and \(\phi =\Vert \varvec{\tau } \Vert\) clearly shows that the stress tensor based indicator leading to \(\alpha (\phi _{t})=\alpha (\Vert \varvec{\tau }_t \Vert )\) exhibits more robust behavior, providing consistent results over wide range of Weissenberg numbers. This kind of time-derivative governed artificial diffusion setup has proved to be quite efficient and safe way of stabilizing the numerical simulations without spoiling the results by excessive (and non-physical) diffusion. One of the drawbacks of this time-derivative driven artificial diffusion stabilization was that when the diffusion coefficient \(\alpha\) was reset at every iteration based on the actual (even very oscillatory) solution, it often led to slow convergence or simulation failure. To eliminate this behavior the floating average of the diffusion coefficient (6) was considered rather than its actual instantaneous value originally defined by (5). This averaging of \(\alpha (\phi _{t})\) over a history of the length L resolved the problem provided that suitable history length L was chosen.

The vanishing diffusion coefficient generalization of the standard artificial stress diffusion technique has proved to be a valuable tool in stabilizing the Oldroyd type model simulations up to moderately high Weissenberg numbers, without nonphysically affecting the final solution. Especially promising is the variant with time-derivative dependent self-adjusting stabilization. The future work in this area should focus on other perspective stabilization variants of this type, with special focus on possible use of such vanishing diffusion stabilization in unsteady simulations.