Flow of viscoelastic fluids past a cylinder at high Weissenberg number: Stabilized simulations using matrix logarithms
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 bywhere p is the pressure, is the velocity vector, the rate-of-deformation tensor is given by
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 . 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 with x in the interval and time such thatwhere a is a constant convection speed and 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 instead of itself. In order to do so, we need an evolution equation for . 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 . The ratio is equal to 2 and the total length of the flow domain is . The flow geometry is shown in Fig. 1. In the following we will use an 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 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 . 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)
- et al.
Constitutive laws for the matrix-logarithm of the conformation tensor
J. Non-Newtonian Fluid Mech.
(2004) - et al.
Simulation of viscoelastic flows using Brownian configuration fields
J. Non-Newtonian Fluid Mech.
(1997) - et al.
A new mixed finite element method for computing viscoelastic flows
J. Non-Newtonian Fluid Mech.
(1995) - et al.
Viscoelastic flow past a confined cylinder of a low density polyethylene melt
J. Non-Newtonian Fluid Mech.
(1997) - et al.
On the discrete EVSS method
Comp. Meth. Appl. Mech. Eng.
(2000) - et al.
A new approach to the deformation fields method for solving complex flows using integral constitutive equations
J. Non-Newtonian Fluid Mech.
(2001) - et al.
Mathematical and physical requirements for successful computations with viscoelastic fluid models
J. Non-Newtonian Fluid Mech.
(1988) Some properties and analytical expressions for plane flow of Leonov and Giesekus models
J. Non-Newtonian Fluid Mech.
(1988)- et al.
Galerkin/least-square finite-element methods for steady viscoelastic flows
J. Non-Newtonian Fluid Mech.
(1999)