Flow of viscoelastic fluids past a cylinder at high Weissenberg number: Stabilized simulations using matrix logarithms

https://doi.org/10.1016/j.jnnfm.2005.01.002Get rights and content

Abstract

The log conformation representation proposed in [R. Fattal, R. Kupferman, Constitutive laws for the matrix-logarithm of the conformation tensor, J. Non-Newtonian Fluid Mech. 123 (2004) 281–285] has been implemented in a FEM context using the DEVSS/DG formulation for viscoelastic fluid flow. We present a stability analysis in 1D and identify the failure of the numerical scheme to balance exponential growth as a possible source for numerical instabilities at high Weissenberg numbers. A different derivation of the log-based evolution equation than in [R. Fattal, R. Kupferman, Constitutive laws for the matrix-logarithm of the conformation tensor, J. Non-Newtonian Fluid Mech. 123 (2004) 281–285] is also presented. We show numerical results for the flow around a cylinder for an Oldroyd-B and a Giesekus model. With the log conformation representation, we are able to obtain solutions beyond the limiting Weissenberg numbers in the standard scheme. In particular, for the Giesekus model the improvement is rather dramatic: there does not seem to be a limit for the chosen model parameter (α=0.01). However, it turns out that although in large parts of the flow the solution converges, we have not been able to obtain convergence in localized regions of the flow. Possible reasons include artefacts of the model and unresolved small scales. More work is necessary, including the use of more refined meshes and/or higher order schemes, before any conclusion can be made on the local convergence problems.

Introduction

The high Weissenberg number problem (HWNP) has been the major stumbling block in computational rheology for the last three decades. The term “HWNP” refers to the empirical observation that all numerical methods break down when the Weissenberg number exceeds a critical value. The precise critical value at which computations break down varies with the problem (including the constitutive model), the method and the mesh used. The precise mechanism that leads to the breakdown has remained something of a mystery (see [2] and the references therein).

A mechanism responsible for at least part of the problems seen at high Weissenberg number has been recently proposed in [3]. In essence, it is a numerical instability caused by the failure to balance the exponential growth of the stress (due to deformation) with convection. Such a failure is common to all methods that approximate the stress by polynomial base functions (the choice of polynomial base function is explicit in finite elements and implicit in finite differences). The remedy proposed in [3] was a change of variables into new variables that scale logarithmically with the stress tensor. Specifically, the constitutive relations were reformulated as equations for the matrix logarithm of the conformation tensor, exploiting the fact that the latter is symmetric positive definite. Numerical experiments for a two-dimensional (2D) lid-driven cavity using a second-order finite difference scheme indicate that schemes based on the new logarithmic formulation are immune to the high Weissenberg number breakdown.

In this paper, we implement the new logarithmic formulation with a finite element method (FEM), and test it for a classical benchmark problem—flow past a cylinder. More specifically, we have in mind the following goals:

  • (1)

    Develop a FEM scheme with much better numerical stability.

  • (2)

    Get more insight into the HWNP by analyzing the numerical breakdown within a FEM point of view.

  • (3)

    Present benchmark results that confirm the validity of the new method at low and moderate Weissenberg numbers, where comparison to existing data is possible.

  • (4)

    Establish new benchmark results at higher Weissenberg numbers.

  • (5)

    Investigate the limitations of the method, and in particular, understand how accuracy is lost at high Weissenberg numbers. Such an understanding is important for the future design of higher order methods.

The paper is structured as follows. In Section 2 we describe the governing equations and introduce the notations used henceforth. In Section 3, we describe the standard FEM discretization based on DEVSS and discontinuous Galerkin. In Section 4, we analyze the stability of the standard discretization, and in particular obtain a simple stability criterion whose violation causes the high Weissenberg number breakdown. In Section 5, we derive the constitutive equation for the matrix logarithm of the conformation tensor; the derivation is based on a different approach than in [1]. In Section 6, we present numerical results for the flow past a cylinder for an Oldroyd-B and a Giesekus model. In particular, we will verify the criterion for breakdown in the standard discretization. We also show the improved stability of the matrix log formulation. A discussion follows in Section 7.

Section snippets

Governing equations

We will consider the flow of viscoelastic fluids at creeping flow conditions (inertia can be neglected) in a spatial region Ω, having a boundary denoted by Γ. For this, we need the following set of equations: the momentum balance, the mass balance for incompressible flows and a constitutive equation describing the stress.

The momentum balance and mass balance are given byp(2ηsD)τ=0,u=0,where p is the pressure, u is the velocity vector, the rate-of-deformation tensor is given by D=12(L+LT

Numerical discretization

For the spatial discretization of the system of equations, we will use the finite element method. We will basically use the same implementation as described in [4] and give a brief summary here.

In order to obtain a proper mixed velocity-stress formulation, we use the Discrete Elastic-Viscous Split Stress (DEVSS) formulation of Guénette and Fortin [5] for the discretization of the linear momentum balance and the continuity equation. For this we introduce an extra variable e=2ηpD. Note that

A stability criterion for exponential profiles

In [3], an analysis of a numerical instability is given for a one-dimensional (1D) problem discretized using a first-order upwind scheme. Here we briefly repeat that and at the same time extend the analysis to our DG scheme.

The ‘toy problem’ we are considering here is as follows: find c(x,t) with x in the interval (0,L) and time t>0 such thatct+acx=bc,where a is a constant convection speed and b>0 is the exponential growth factor. Eq. (14) can be interpreted as a 1D representation of the

Evolution equation for the logarithm of the conformation tensor

In the previous section, we have identified the failure to resolve exponential profiles as a major restriction on the stability of standard schemes. It is suggested that this restriction can be lifted by solving for the matrix-logarithm of the conformation tensor c instead of c itself. In order to do so, we need an evolution equation for s=logc. In [1] this evolution equation has been derived by a decomposition of the velocity gradient. In this paper, we give an alternative derivation using an

Problem description

We consider the planar flow past a cylinder of radius R positioned between two flat plates separated by a distance 2H. The ratio H/R is equal to 2 and the total length of the flow domain is 30R. The flow geometry is shown in Fig. 1. In the following we will use an (x,y) co-ordinate system with the origin positioned at the center of the cylinder.

We assume the flow to be periodic. This means that we periodically extend the flow domain such that cylinders are positioned 30R apart. This avoids

Conclusions and discussion

It has been shown that also in the FEM implementation (DEVSS with Discontinuous Galerkin), the log conformation representation removes the catastrophic breakdown present in the standard FEM implementation. We use a standard benchmark problem: the flow around a cylinder using an Oldroyd-B model and also tested a Giesekus model with α=0.01. Especially the Giesekus model shows a dramatic improvement of the numerical stability. That doesn’t mean all problems are solved.

It turns out that high

Acknowledgements

We thank Frank Baaijens for stimulating discussions and for giving the first author the opportunity to present this work at the ICR2004.

This research was funded in part by the Applied Mathematical Sciences subprogram of the Office of Energy Research of the US Department of Energy under Contract DE-AC03-76-SF00098.

References (18)

There are more references available in the full text version of this article.

Cited by (0)

View full text